arXiv:1505.03786vl [astro-ph.GA] 14 May 2015
Quasar Quartet Embedded in Giant Nebula
Reveals Rare Massive Structure in Distant Universe
Joseph F. Hennawi 1 *, J. Xavier Prochaska 2 ,Sebastiano Cantalupo 2,3 ,
Fabrizio Arrigoni-Battaia 1
1 Max-Planck-Institut fur Astronomie, Konigstuhl, Germany
2 University of California Observatories-Lick Observatory, UC Santa Cruz, CA
'ETH Zurich, Institute of Astronomy, Zurich, Switzerland
*To whom correspondence should be addressed; E-mail: firstname.lastname@example.org
All galaxies once passed through a hyper-luminous quasar phase powered by
accretion onto a supermassive black hole. But because these episodes are brief,
quasars are rare objects typically separated by cosmological distances. In a
survey for Lyman-o emission at redshift z & 2, we discovered a physical as¬
sociation of four quasars embedded in a giant nebula. Located within a sub¬
stantial overdensity of galaxies, this system is likely the progenitor of a massive
galaxy cluster. The chance probability of finding a quadruple-quasar is esti¬
mated to be ~ Hr 7 , implying a physical connection between Lyman-a nebulae
and the locations of rare protoclusters. Our findings imply that the most mas¬
sive structures in the distant Universe have a tremendous supply (~ 10 n solar
masses) of cool dense (volume density ~ 1 cm -3 ) gas, in conflict with current
Cosmologists do not fully understand the origin of supermassive black holes (SMBHs) at the
centers of galaxies, and how they relate to the evolution of the underlying dark matter, which
forms the backbone for structure in the Universe. In the current paradigm, SMBHs grew in
every massive galaxy during a luminous quasar phase, making distant quasars the progenitors
of the dormant SMBHs found at the centers of nearby galaxies. Tight correlations between the
masses of these local SMBHs and both their host galaxy ( 1 ) and dark matter halo masses (2)
support this picture, further suggesting that the most luminous quasars at high-redshift should
reside in the most massive galaxies. It has also been proposed that quasar activity is triggered
by the frequent mergers that are a generic consequence of hierarchical structure formation (3,4).
Indeed, an excess in the number of small-separation binary quasars (5, 6), as well as the mere
existence of a handful of quasar triples (7,8), support this hypothesis. If quasars are triggered by
mergers, then they should preferentially occur in rare peaks of the density field, where massive
galaxies are abundant and the frequency of mergers is highest (9).
Following these arguments, one might expect that at the peak of their activity z ~ 2 — 3,
quasars should act as signposts for protoclusters, the progenitors of local galaxy clusters and
the most massive structures at that epoch. However, quasar clustering measurements (10,11)
indicate that quasar environments at z ~ 2 — 3 are not extreme: these quasars are hosted by
dark matter halos with masses M halo ~ 10 12 5 M 0 (where M e is the mass of the Sun), too
small to be the progenitors of local clusters (12). But the relationship between quasar activity
and protoclusters remains unclear, owing to the extreme challenge of identifying the latter at
high-redshift. Indeed, the total comoving volume of even the largest surveys for distant galaxies
at z ~ 2 — 3 is only ~ 10' Mpc 3 , which would barely contain a single rich cluster locally.
Protoclusters have been discovered around a rare-population of active galactic nuclei (AGN)
powering large-scale radio jets, known as high-redshift radio galaxies (HzRGs) (13). The
HzRGs routinely exhibit giant ~ 100 kpc nebulae of luminous Lyo emission L Ly a ~ 10 44 ergs 4 .
Nebulae of comparable size and luminosity have also been observed in a distinct population of
objects known as ‘Lycc blobs’ (LABs) (14,15). The LABs are also frequently associated with
AGN activity (16-18), although lacking powerful radio jets, and appear to reside in overdense,
protocluster environments (15,19, 20). Thus among the handful of protoclusters (21) known,
most appear to share two common characteristics: the presence of an active SMBH and a giant
We have recently completed a spectroscopic search for extended Lycc emission around a
sample of 29 luminous quasars at z ~ 2, each the foreground (f/g) member of a projected
quasar pair (22). Analysis of spectra from the background (b/g) members in such pairs reveals
that quasars exhibit frequent absorption from a cool, metal-enriched, and often optically thick
medium on scales of tens to hundreds kpc (22-28). The UV radiation emitted by the luminous
f/g quasar can, like a flashlight, illuminate this nearby neutral hydrogen, and power a large-scale
Lya-cmission nebula, motivating our search. By construction, our survey selects for exactly
the two criteria which seem to strongly correlate with protoclusters: an active SMBH and the
presence of a large-scale Lyo emission nebula.
Of the 29 quasars surveyed, only SDSSJ0841+3921 exhibited extended large-scale (> 50 kpc)
Lycc emission above our characteristic sensitivity limit of 6 x 10~ 18 erg s -1 cm -2 arcsec -2 (2a).
We designed a custom narrow-band filter tuned to the wavelength of Lyo at the f/g quasar red-
shift z = 2.0412 (A C e n ter = 3700A, FWHM A = 33A), and imaged the field with the Keck/LRIS
imaging spectrometer for 3hrs on UT 12 November 2012. The combined and processed images
reveal Lya emission from a giant filamentary structure centered on the f/g quasar and extending
continuously towards the b/g quasar (see Fig. 1). This nebulosity has an end-to-end size of 37"
corresponding to 310 kpc and a total line luminosity ^Lya = 2.1 x 10 44 ergs _1 , making it one
of the largest and brightest nebula ever discovered.
The giant nebula is only one of the exceptional properties of SDSSJ 08414-3921. Our
images reveal three relatively compact candidate Lyct emitting sources with faint continuum
magnitudes V ~ 23 — 24, embedded in the Lyct filament and roughly aligned with its major
axis. Follow-up spectroscopy reveals that the sources labeled AGN1, AGN2, and AGN3 are
three AGN at the same redshift as the f/g quasar (see right panel of Fig. land (29)), making this
system the only quadruple AGN known. Adopting recent measurements of small-scale quasar
clustering (30), we estimate that the probability of finding three AGN around a quasar with
such small separations is ~ 10~ 7 (31). Why then did we discover this rare coincidence of AGN
in a survey of just 29 quasars? Did the giant nebula mark the location of a protocluster with
dramatically enhanced AGN activity?
To test this hypothesis, we constructed a catalog of Lyo;-emitting galaxies (LAEs), and
computed the cumulative overdensity profile of LAEs around SDSSJ0841+3921, relative to
the background number expected based on the LAE luminosity function (32) (see Fig. 2). To
perform a quantitative comparison to other giant Lya nebulae, many of which are known to
coincide with protoclusters, we measured the giant nebulae-LAE cross-correlation function for
a sample of eight systems - six HzRGs and two LABs - for which published data was available
in the literature (33). In Fig. 2, we compare the overdensity profile around SDSSJ0841+3921 to
this giant nebulae-LAE correlation function. On average, the environment of HzRGs and LABs
hosting giant Lya nebulae (red line) is much richer than that of radio-quiet quasars (10) (blue
line), confirming that they indeed reside in protoclusters. Furthermore, the clustering of LAEs
around SDSSJ0841+3921 has a steeper overdensity profile, and exceeds the average protoclus¬
ter by a factor of > 20 for R < 200 kpc and by ~ 3 on scales of R ~ 1 Mpc. In addition to
the overdensity of four AGN, the high number of LAEs surrounding SDSSJ0841+3921 make
it one of the most overdense protoclusters known at z ~ 2 — 3.
The combined presence of several bright AGN, the Lyo: emission nebula, and the b/g absorp¬
tion spectrum, provide an unprecedented opportunity to study the morphology and kinematics
of the protocluster via multiple tracers, and we find evidence for extreme motions (34). Specif¬
ically, AGN 1 is offset from the precisely determined systemic redshift (35) of the f/g quasar by
+1300 ± 400 km s’ 1 . This large velocity offset cannot be explained by Hubble expansion - the
miniscule probability of finding a quadruple quasar in the absence of clustering P rsj 10~ 13 (31)
and the physical association between the AGN and giant nebula demand that the four AGN
reside in a real collapsed structure - and thus provides an unambiguous evidence for extreme
gravitational motions. In addition, our slit spectra of the giant Lya nebula reveal extreme kine¬
matics of diffuse gas (Fig. 3), extending over a velocity range of —800 to +2500 km s -1 from
systemic. Furthermore, there is no evidence for double-peaked velocity profiles characteristic
of resonantly-trapped Lya, which could generate large velocity widths in the absence of corre¬
spondingly large gas motions. Absorption line kinematics of the metal-enriched gas, measured
from the b/g quasar spectrum at impact parameter of R± = 176 kpc (Fig. 4), show strong ab¬
sorption at « +650 km s^ 1 with a significant tail to velocities as large as ~ 1000 km s -1 . It is
of course possible that the extreme gas kinematics, traced by Lyo emission and metal-line ab¬
sorption, are not gravitational but rather arise from violent large-scale outflows powered by the
multiple AGN. While we cannot completely rule out this possibility, the large velocity offset of
+1300 km s’ 1 between the f/g quasar and the emission redshift of AGN1 can only result from
One can only speculate about the origin of the dramatic enhancement of AGN in the
SDSSJ0841+3921 protocluster. Perhaps the duty cycle for AGN activity is much longer in pro¬
toclusters, because of frequent dissipative interactions (5,6), or an abundant supply of cold gas.
A much larger number of massive galaxies could also be the culprit, as AGN are known to trace
massive halos at z ~ 2. Although SDSSJ0841+3921 is the only example of a quadruple AGN
with such small separations, previously studied protoclusters around HzRGs and LABs, also
occasionally harbor multiple AGN (13,36). Regardless, our discovery of a quadruple AGN and
protocluster from a sample of only 29 quasars suggests a link between giant Lya nebulae, AGN
activity, and protoclusters - similar to past work on HzRGs and LABs - with the exception that
SDSSJ 0841+3921 was selected from a sample of normal radio-quiet quasars. From our survey
and other work (37), we estimate that ~ 10% of quasars exhibit comparable giant Lya nebulae.
Although clustering measurements imply that the majority of z ~ 2 quasars reside in moderate
overdensities (10-12), we speculate that this same 10% trace much more massive protoclus¬
ters. SDSSJ0841+3921 clearly supports this hypothesis, as does another quasar-protocluster
association (10,38), around which a giant Lya nebula was recently discovered (39,40).
Given our current theoretical picture of galaxy formation in massive halos, an association
between giant Lya nebulae and protoclusters is completely unexpected. The large Lya lumi¬
nosities of these nebulae imply a substantial mass ( rs./ 10 11 M 0 ) of cool (T ~ 10 4 K) gas (41),
whereas cosmological hydrodynamical simulations indicate that already by z ~ 2 — 3, baryons
in the massive progenitors (M ha i 0 > 10 13 M 0 ) of present-day clusters are dominated by a hot
shock-heated plasma T ~ 10' K (42,43). These hot halos are believed to evolve into the X-ray
emitting intra-cluster medium observed in clusters, for which absorption-line studies indicate
negligible < 1% cool gas fractions (44). Clues about the nature of this apparent discrepancy
come from our absorption line studies of the massive ~ 10 12 ' 5 M e halos hosting z ~ 2 — 3
quasars. This work reveals substantial reservoirs of cool gas > 10 10 M e (22-28), manifest as
a high covering factor ~ 50% of optically thick absorption, several times larger than predicted
by hydrodynamical simulations (42,43). This conflict most likely indicates that current simula¬
tions fail to capture essential aspects of the hydrodynamics in massive halos at z ~ 2 (27,42),
perhaps failing to resolve the formation of clumpy structure in cool gas (41).
If illuminated by the quasar, these large cool gas reservoirs in the quasar circumgalactic
medium (CGM) will emit fluorescent Lya photons, and we argue that this effect powers the neb¬
ula in SDSSJ0841+3921 (45). But according to this logic, nearly every quasar in the Universe
should be surrounded by a giant Lya nebulae with size comparable to its CGM (~ 200 kpc).
Why then are these giant nebulae not routinely observed?
This apparent contradiction can be resolved as follows. If this cool CGM gas is illumi¬
nated and highly ionized, it will fluoresce in the Lya line with a surface brightness scaling
as SB LyQ oc A H nn, where Ah is the column density of cool gas clouds which populate the
quasar halo, and % is the number density of hydrogen atoms inside these clouds. Note the total
cool gas mass of the halos scales as M cool oc f? 2 A H , where R is the radius of the halo (45).
Given our best estimate for the properties of the CGM around typical quasars (% ~ 0.01 cm -3
and A'n ~ 10 2 ° cm 2 or M coo \ ~ 10 11 M 0 ) (22, 26, 27), we expect these nebulae to be ex¬
tremely faint SB Lya ~ 10~ 19 ergs -1 cm -2 arcsec -2 , and beyond the sensitivity of current in¬
struments (22). One comes to a similar conclusion based on a full radiative transfer calculation
through a simulated dark-matter halo with mass M hak) « 10 12,5 M 0 (41). Thus the factor of
~ 100 times larger surface brightness observed in the SDSSJ0841+3921 and other protocluster
nebulae, arises from either a higher %, iV H (and hence higher M coo 0, or both. The cool gas
properties required to produce the SDSSJ0841+3921 nebula can be directly compared to those
deduced from an absorption line analysis of the b/g quasar spectrum (46).
The b/g quasar sightline pierces through SDSSJ0841+3921 at an impact parameter of R± =
176 kpc, giving rise to the absorption spectrum shown in Fig. 4. Photoionization modeling of
these data constrains the total hydrogen column density to be \og l0 Nn = 20.4 ± 0.4 (45),
implying a substantial mass of cool gas 1.0 x 10 11 M 0 < M cooi < 6.5 x 10 11 M 0 within
r = 250 kpc. Assuming that the Ly a emitting gas has the same column density as the gas
absorbing the b/g sightline, reproducing the large fluorescent Lya surface brightness requires
that this gas be distributed in compact r c i ou d ~ 40 pc clouds at densities characteristic of the
interstellar medium n\\ ~ 2 cm' 3 , but on ~ 100 kpc scales in the protocluster.
Clues to the origin of these dense clumps of cool gas comes from their high enrichment level,
which we have determined from our absorption line analysis (46) to be greater than l/10th of
the solar value. At first glance, this suggests that strong tidal interactions due to merger activity
or outflows due to powerful AGN feedback are responsible for dispersing dense cool gas in the
protocluster. However, the large cool gas mass ~ 10 11 M 0 and high velocities ~ 1000 km s -1 ,
imply an extremely large kinetic luminosity L win( j ~ 10 44 6 for an AGN powered wind, making
the feedback scenario implausible (25). An even more compelling argument against a merger or
feedback origin comes from the extremely small cloud sizes r cloud ~ 40 pc implied by our mea¬
surements. Such small clouds moving supersonically ~ 1000 km s _1 through the hot T ~ 10' K
shock-heated plasma predicted to permeate the protocluster, will be disrupted by hydrodynamic
instabilities in ~ 5 x 10 6 yr, and can thus only be transported ~ 5 kpc (47). These short dis¬
ruption timescales instead favor a scenario where cool dense clouds are formed in situ, perhaps
via cooling and fragmentation instabilities, but are short-lived. The higher gas densities might
naturally arise if hot plasma in the incipient intra-cluster medium pressure confines the clouds,
compressing them to high densities (48, 49). Emission line nebulae from cool dense gas has
also been observed at the centers of present-day cooling flow clusters (50, 51), albeit on much
smaller scales < 50 kpc. The giant Lya nebulae in z ~ 2 — 3 protoclusters might be manifesta¬
tions of the same phenomenon, but with much larger sizes and luminosities, reflecting different
physical conditions at high-redshift.
The large reservoir of cool dense gas in the protocluster SDSSJ0841+3921, as well as
those implied by the giant nebulae in other protoclusters, appear to be at odds with our current
theoretical picture of how clusters form. This is likely symptomatic of the same problem of too
much cool gas in massive halos already highlighted for the quasar CGM (27,42,43). Progress
will require more cosmological simulations of massive halos M > 10 13 M 0 at z ~ 2, as well
as idealized higher-resolution studies. In parallel, a survey for extended Lya emission around
~ 100 quasars would uncover a sample of ~ 10 giant Lya nebulae likely coincident with
protoclusters, possibly also hosting multiple AGN, and enabling continued exploration of the
relationship between AGN, cool gas, and cluster progenitors.
Fig. 1: Narrow and broad band images of the field surrounding SDSSJ0841+3921 Left:
Continuum-subtracted, narrow-band image of the field around f/g quasar. The color map and the contours
indicate the Lya surface brightness (upper color bar) and the signal-to-noise ratio per arcsec 2 aperture
(lower color bar), respectively. This image reveals a giant Lya nebula on the northern side of the f/g
quasar and several compact, bright Lya emitters in addition to the f/g quasar. Three of these have been
spectroscopically confirmed as active galactic nuclei (AGN) at the same redshift. Right: Corresponding
V -band continuum image of the field presented at left with the locations of the four AGN marked. The
AGN are roughly oriented along a line coincident with the projected orientation of the Lya nebula. We
also mark the position of the b/g quasar, which is not physically associated with the quadruple AGN
system, but whose absorption spectrum probes the gaseous environment of the f/g quadruple AGN and
protocluster (see Fig. 4).
Fig. 2: Characterization of the protocluster environment around SDSSJ0841+3921 . The
data points indicate the cumulative overdensity profile of LAEs <5(< R) as a function of impact parameter
R from the f/g quasar in SDSSJ0841+3921, with Poisson error bars. The red curve shows the predicted
overdensity profile, based on our measurement of the giant nebulae-LAE cross-correlation function de¬
termined from a sample of eight systems - six HzRGs and two LABs - for which published data was
available in the literature. Assuming a power-law form for the cross-correlation £ cross = (r/?’o) -7 , we
measure the correlation length ro = 29.3 ± 4.9 h ~ 1 Mpc, for a fixed value of 7 = 1.5. The gray shaded
region indicates the 1 a error on our measurements based on a bootstrap analysis, where both ro and
7 are allowed to vary. The solid blue line indicates the overdensity of Lyman break galaxies (LBGs)
around radio-quiet quasars based on recent measurements (10), with the dotted blue lines the ler error
on this measurement. On average, the environment of HzRGs and LABs hosting giant Lya nebulae is
much richer than that of radio-quiet quasars (10), confirming that they indeed reside in protoclusters.
SDSSJ0841+3921 exhibits a dramatic excess of LAEs compared to the expected overdensities around
radio-quiet quasars (blue curve). Its overdensity even exceeds the average protocluster (red curve) by a
factor of > 20 for R < 200 kpc decreasing to an excess of ~ 3 on scales of II ~ 1 Mpc, and exhibits a
much steeper profile.
Fig. 3: Ly a spectroscopy of the giant nebula and its associated AGN. Upper Left: The
spectroscopic slit locations (white rectangles) for three different slit orientations are overlayed on the
narrow band image of the giant nebula. The locations of the f/g quasar (brightest source), b/g quasar,
AGN1, and AGN2 are also indicated. Two-dimensional spectra for Slitl (Upper Right), Slit2 (Lower
Left), and Slit3 (Lower Right) arc shown in the accompanying panels. In the upper right and lower
left panels, spatial coordinates refer to the relative offset along the slit with respect to the f/g quasar.
Spectra of AGN1 are present both in Slit 1 (upper right) and Slit 3 (lower right) at spatial offsets 75kpc
and 25kpc, respectively, while the Lya spectrum of AGN2 is located at a spatial offset — 60kpc in Slit 1
(upper right). The b/g quasar spectra are located in both Slit 2 (lower left) and Slit 3 (lower right) at the
same spatial offset of 176 kpc. The spectroscopic observations demonstrate the extreme kinematics of
the system: AGN 1 has a velocity of +1300 ± 400 km s _1 relative to the f/g quasar and the Lycc emission
in the nebula exhibits motions ranging from —800km s” 1 (at ~ 100 kpc offset in Slit 3, lower right) to
+2500km 1 (at ~ 100 kpc in Slit 1, upper right). A 3 X 3 pixel boxcar smoothing, which corresponds
120 km s -1 X 0.8", has been applied to the spectra. In each two-dimensional spectrum, the zero velocity
corresponds to the systemic redshift of the f/g quasar. The color bars indicate the flux levels in units of
erg s -1 cm -2 arcsec -2 A -1 . 9
Fig. 4: Absorption line spectrum of cool gas in SDSSJ0841+3921 . Spectrum of the absorbing
gas detected in the b/g quasar sightline at an impact parameter of 176 kpc from the f/g quasar. The
gas shows strong Hi and low-ionization state metal absoiption, offset by ~ 650km s -1 from the f/g
quasar’s systemic redshift. The CII absoiption in particular exhibits a significant tail to velocities as
large as ~ 1000 km s -1 , providing evidence for extreme gas kinematics. We modeled the strong HI
Lya absorption with a Voigt profile (blue curve with grey band indicating uncertainty) and estimate a
column density log Api = 19.2 ± 0.3. The strong low and intermediate ion absorption (Sill, CII, Silll)
and correspondingly weak high-ion absorption (CIV, SilV) indicate that the gas is not highly ionized,
and our photoionization modeling (45) implies log 10 rrni = —1.2 ± 0.3 or log 10 Ah = 20.4 ± 0.4. We
estimate a conservative lower-limit on the gas metallicity to be 1/10 of the solar value.
We thank the staff of the W.M. Keck Observatory for their support during the installation and
testing of our custom-built narrow-band filter. We are grateful to B. Venemans and M. Prescott
for providing us with catalogs of LAE positions around giant nebulae in electronic format.
We also thank the members of the ENIGMA group (http://www.mpia-hd.mpg.de/ENIGMA/)
at the Max Planck Institute for Astronomy (MPIA) for helpful discussions. JFH acknowledges
generous support from the Alexander von Humboldt foundation in the context of the Sofja Ko-
valevskaja Award. The Humboldt foundation is funded by the German Federal Ministry for
Education and Research. J.X.P. acknowledge support from the National Science Foundation
(NSF) grant AST-1010004. The data presented here were obtained at the W.M. Keck Observa¬
tory, which is operated as a scientific partnership among the California Institute of Technology,
the University of California and NASA. The Observatory was made possible by the financial
support of the W.M. Keck Foundation. We acknowledge the cultural role that the summit of
Mauna Kea has within the indigenous Hawaiian community. We are most fortunate to have the
opportunity to conduct observations from this mountain. The data reported in this paper are
available through the Keck Observatory Archive (KOA).
Supplementary Online Material
Materials and Methods
Figs. SI to S10
Tables S1 to S6
Supplementary Online Material
This PDF file includes
Materials and Methods
Figs. SI to S10
Tables S1 to S6
1 Optical Observations
1.1 Discovery Spectra
The source SDSSJ084158.47+392121.0 (f/g quasar) was targeted as a quasar by the Sloan Digi¬
tal Sky Survey (SDSS) through their color-selection algorithms and observed spectroscopically
in the standard survey (52). It has a cataloged emission redshift of z — 2.046. Our analysis of
the SDSS imaging data revealed a second, neighboring source at SDSSJ084159.26+392140.0
(b/g quasar) with colors characteristic of a z ~ 2 quasar, characterizing this system as a can¬
didate quasar pair. In the course of our ongoing spectroscopic campaign to discover close
quasar pairs at z > 2 (6, 53), we confirmed this source to have z = 2.214 implying a pro¬
jected quasar pair with a physical separation of R± = 176 kpc at the f/g quasar redshift.
SDSSJ084159.26+392140.0 was also observed by the SDSS-III survey in the Baryonic Oscil¬
lating Spectroscopic Survey campaign at a spectral resolution of R & 2000 and with wavelength
coverage A ~ 3600 — 10, 000A (54).
On UT 2007 Jan 18, we observed the quasar pair with the Low Resolution Imaging Spec¬
trograph (LRIS) (55). These data were taken to study the intergalactic medium probed by the
quasar pair, to examine the HI gas associated with the circumgalactic medium of the f/g quasar,
and to search for fluorescent emission associated with the f/g quasar. The latter analyses of
these data have been presented in previous works (22, 26-28). Summarizing the observations,
we used LRIS in multi-slit mode with a custom designed slitmask which allowed placement
of one slit on the known quasars and other slits on additional quasar candidates in the field.
Specifically, a slit was placed at the position angle PA=25.8° between the f/g and b/g quasars,
allowing them to be observed simultaneously. LRIS is a double spectrograph with two arms
giving simultaneous coverage of the near-UV and optical. We used the D460 dichroic with the
1200 lines mm -1 grism blazed at 3400 A on the blue side, resulting in wavelength coverage
of « 3250 — 4300 A. The dispersion of this grism is 0.50 A per pixel and our 1" slits give a
resolution with full width at half maximum (FWHM) FWHM~ 160km s~ 1 . On the red side we
used the R1200/5000 (covering ~ 4700 —6000A) and R300/5000 (covering ~ 4700 — 10, 000A)
gratings having a FWHM of « 100km s -1 and 400km s’ 1 respectively.
The science frames were complemented by a series of calibration images: arc lamp, dome
flat, twilight sky, and standard star spectra with the same instrument configurations. All of these
exposures were reduced using the LowRedux (http://www.ucolick.org/~xavier/LowRedux/)
pipeline which bias subtracts and flat fields the images, corrects for non-uniform illumination,
derives a wavelength solution, performs sky subtraction, optimally extracts the sources, and
fluxes the resultant spectra. The ID spectra are corrected for instrument flexure and shifted to a
heliocentric, vacuum-corrected system.
Using custom software, we also coadded the 2D spectral images to search for diffuse Lyct
emission surrounding the f/g quasar (22). This software enables us to model and subtract the
spectral PSF of sources in the 2D spectra. Extended Lycc emission will be manifest as residual
flux in our 2D sky-and-PSF-subtracted images which is inconsistent with being noise. To vi¬
sually assess the statistical significance of any putative emission feature, we define a x image
Xsky+PSF = (DATA — SKY — OBJECTS) /a. If our model is an accurate description of the
data, the distribution of pixel values in the y s k y +PSF should be a Gaussian with unit variance.
The middle row of the images in Fig. SI shows this quantity for SDSSJ 0841+3921 for each
slit orientation. The lower row of these images shows y s k y = (DATA — SKY) /a. The upper
row shows a smoothed version of the middle row, helpful for identifying extended emission.
Specifically, the smoothed images are given by yy mth = CONx - objects] , w ^ ere
the CONVOL operation denotes smoothing of the stacked images with a symmetric Gaus¬
sian kernel (same spatial and spectral widths) with FWHM=235 km s -1 (dispersion cr smt h =
100 km s -1 ), which is 5.7 pixels, or 1.4 times the spectral resolution element, and corresponds
to FWHM = 1.5" spatially. The operation CONVOL 2 represents convolution with the square
of the smoothing kernel. With this definition of yy mtll the distribution of the pixel values in the
smoothed image will still obey Gaussian statistics, although they are of course correlated, and
hence not independent.
Sky and object PSF-subtracted y-maps for all of the slit-orientations that we used to char¬
acterize extended emission in SDSSJ 0841+3921 are shown in Fig. SI. Fig. 3 of the Main text,
shows just the smoothed sky-subtracted images with a color-map chosen to accentuate the faint
extended emission and a slightly different smoothing, namely a 3 x 3 pixel boxcar smoothing,
which corresponds 120 km s 1 x 0. 8". The y-maps in Fig. SI enable the reader to objectively
assess the statistical significance of all emission features in the unsmoothed data. Note that in
all of these maps, only the PSF model of the f/g and b/g quasars have been subtracted, whereas
neither AGN1 nor AGN2 were removed from the other slit-orientations (see Fig. 1 in main text).
Following the calibration procedure described in (22), we deduce the following spectroscopic
surface brightness limits (la) for the Fya line at the f/g quasar redshift z = 2.0412: Slitl,
SBi ct = 2.2 x 10 -18 ergs -1 cm -2 arcsec -2 ; Slit2, SBi ct = 2.3 x 10 -18 ergs -1 cm -2 arcsec -2 ;
and Slit3, SB i a = 4.8 x 10 18 ergs 1 cm 2 arcsec 2 , where these SBs are computed in win¬
dows of 700 km s -1 x 1.0", which corresponds to an aperture of 700 km s -1 x 1.0 arcsec 2 on the
sky because we always used a 1.0" slit. The depth that we attain in our Fya spectroscopy is com¬
parable to that achieved by our narrow-band imaging SBi ct = 1.7 x 10 -18 ergs -1 cm -2 arcsec -5
(see next section).
The original observation from the Quasars Probing Quasars (QPQ) (22) survey corresponds
to slit-orientation Slit2 (see Fig. SI and Fig. 3 in main text). Our initial visual inspection
of the Lya emission map revealed a bridge of Lya emission along the slit connecting the
quasar pair. Although the SB varies with position, this structure has a characteristic SB Lya ~
10 -1 ' ergs -1 cm -2 arcsec -2 and is detected at ~ 4 — 5a at location R± ~ +100 kpc between
the two quasars. The emission extends along the slit from —80 kpc below the f/g quasar to
+150 kpc between the two quasars, and possibly to +250 kpc. With an end-to-end size of
~ 250 kpc, this represented one of the largest Lyo: emission nebula ever detected, and thus
motivated the additional imaging and spectroscopic observations, which are the focus of this
1.2 Narrow Band Imaging
We purchased a custom-designed narrow-band (NB) filter from Andover Corporation, sized
to fit within the grism holder of the Keck/LRISb camera. The filter was tuned to A cen ter =
3700A, and designed with a narrow band-pass FWIIM A = 33A to minimize sky background
while maintaining throughput. On UT 12 November 2012, we imaged the ~ 5' x 7' field-of-
view surrounding SDSSJ0841+3921, offset to place the quasar pair on the CCD with highest
quantum efficiency. We observed for a total of 3.2 hours in a series of dithered, 1280s exposures.
Conditions were clear with sub-arcsecond atmospheric seeing. In parallel, we obtained 3hrs of
broad-band V images with the LRISr camera. The instrument was configured with the D460
dichroic and the detectors of the blue camera were binned 2x2 to minimize read noise.
The images were reduced using standard routines within the IRAF reduction software pack¬
age. This includes bias subtraction, flat fielding and an illumination correction. A combination
of twilight sky flats and unregistered science frames were used to produce flat-field images and
illumination corrections in each band. The individual frames were sequentially registered to the
SDSS-DR7 catalog using the SExtractor ( 56) and SCAMP (57) packages. The RMS uncertainty
in the astrometry of our registered images is approximately 0.2". Finally, the corrected frames
from each band were average-combined using the SWarp package (55).
We calibrated the photometry of our images as follows. At the beginning and end of the
night, we observed two spectrophotometric stars (G191b2b and Feige34) with the NB filter
under clear conditions. Neither star has significant spectral features in the relevant wavelength
range covered by our NB filter. For the broad-band images, we observed the standard star
field PG0231+051. To compute the zero-point for the narrow-band images, we compared the
measured count rates of Feige34 and G191b2b with the expected fluxes estimated by convolving
the standard star spectra (resolution of lA) (59) with the normalized filter transmission curve.
We calculate an average, zero-point magnitude of 23.58 mag, with a few percent difference
between the two stars. With this calibration, we deduce that the la limiting surface brightness
of our combined NB images is SBi ct = 1.7 x 10 -18 ergs -1 cm -2 arcsec -2 for a 1.0 arcsec 2
For the broad-band images, we compared the number of counts per second of the five stars in
the PG0231+051 -field with their tabulated R-band magnitudes (60). The derived zero-point for
the five stars are consistent to within a few percent and we adopt the average value: V Z p = 28.07.
Because the Feige34 and the PG0231+051 fields were observed with an airmass (AMps 1.2)
similar to our science field, we did not correct the individual images before combination. We
estimate that the correction would be on the order of few percent. The la limiting point source
magnitude for our combined broad band images is V = XX.
To isolate the emission in the Lya line we estimated and then subtracted the continuum
emission underlying the NB3700 filter. We estimate the continuum using the 17-band (not
affected by the Lya line) and assuming a flat continuum slope (in frequency), i.e. (3\ = —2.
This assumption is dictated by our dataset, i.e. we can rely only on a deep 17-band image to
estimate the continuum, and it follows the prescriptions of previous work (61) which showed
that /3\ ~ —1.9 for the most luminous sources and higher values for the faintest. Thus, assuming
a flat continuum slope gives us a conservative estimate of the Lya emission and its equivalent
width. As the V --band image has slightly worse seeing (~ 1.3"), we convolved the NB image to
the same seeing before subtracting the continuum. The continuum subtraction has been applied
using the following formula
FWHM N B3700 TrNB3700 y
FWHMy Try ’
( 1 )
where Ly a is the final subtracted image, NB3700 is the narrow-band image, 17 is the 17-band
image, and Tr NB 37 oo an d Try are the peak transmission values for the NB3700 and 17-band
filters, respectively. The result of this procedure is a Lya only narrow-band image, which is
shown in the left panel of Fig. 1.
1.3 Follow-up Spectroscopy
After the discovery of an extended nebula in our NB imaging, we opted to obtain additional
spectroscopy of both the diffuse nebula and several compact sources near the projected quasar
pair during the same observing run. On UT 13 November 2012 in clear observing conditions
and seeing varying from FWHM = 0.6 — 1.0", we used the LRIS spectrometer with a 1.0"
longslit, configured with the D460 dichroic, B1200/3400 grism, and R600/7500 grating. We
obtained two additional orientations of the longslit as illustrated in Fig. 3: (Slit 1) along the line
connecting source f/g quasar and AGN2 with a PA=227°; (Slit 3) along the line connecting the
b/g quasar and AGN1 with a PA=343°. We exposed for a total of 4800s and 2400s for Slit 1
and Slit 3, respectively. The PA of Slit 1 was chosen to be aligned with the f/g quasar and
a nearby radio source detected in the FIRST survey (see Supp. 2.2), which we learned about
just prior to the November 2012 Keck observations. This turned out to be an AGN at the same
redshift (AGN2). The slit orientation also conveniently allowed us to simultaneously observe a
bright spot in the Lya image, which is also the companion AGN 1. The orientation of Slit 3 was
chosen to be aligned with the b/g quasar as a position reference and two other bright spots in
the nebula, one of which was also covered by Slit 1, i.e. AGN1. Coordinates, photometry, and
other information about these targets are listed in Table S1. Corresponding calibration frames
were obtained and these data were reduced with the same procedures described above. The slit
orientation of our original January 2007 spectroscopy is also indicated in Fig. 3 as Slit 2, which
simultaneously observed the f/g and b/g quasars at PA=25.8° for an exposure time of 1800s (but
in better conditions).
Our Lya image (left panel of Fig. 1) reveals several compact Lyo-cmitting sources (LAEs)
with rest-frame Wj jya > 20A in the vicinity of the f/g quasar, which are hence very likely to be at
the same redshift. On UT 2012 Dec 14, we targeted two of these LAEs with the DEIMOS spec¬
trometer (62) on the Keck II telescope in partly cloudy conditions. Specifically, we employed
the 1.0" long-slit mask oriented to cover the object labeled AGN3 in Fig. 1 (see also Table SI),
and another LAE which we designate as Targetl at RA=08:41:58.8, DEC= +39:21:57.4,
which lies off of the image in Fig. 3. The instrument was configured with the 600ZD grat¬
ing tilted to a central wavelength of 7200A which provides coverage from « 4600 — 9800A
with a spectral resolution FWFIM^ 235km s -1 . We took two exposures totaling 4500s. These
data were reduced and extracted with the SPEC2D pipeline (63, 64) and fluxed using a spec¬
trum of the spectrophotometric standard star Feige 34 taken that (non-photometric) night with
the same instrument configuration. Note that given the limited blue sensitivity of DEIMOS,
these spectra did not cover the Lya emission line at ~ 3700A. Owing to its faint continuum
magnitude (V = 25.2 mag) the spectrum of Targetl was inconclusive, although its large Lya
equivalent width I+r jyQ = 97A suggests it is a real LAE. A setup covering Lya is likely neces¬
sary to spectroscopically confirm this source. The object AGN3 is an AGN at the same redshift
as the f/g quasar, as described in the next section.
2 Discovery of Four AGN
Our follow-up spectroscopy reveals three additional AGN within 6 < 18" or R± < 150 kpc
of the f/g quasar, with very similar redshifts. These objects are labeled as AGN1-3 in Fig. 3,
and their optical spectra are shown in Fig. S2. Relevant information for all five of the AGN
associated with the nebula, namely the f/g quasar, the three AGN discovered at the same redshift,
and the b/g quasar are provided in Table S1.
Motivated by the discovery of these AGN, we force-photometered the SDSS and WISE
survey images at the coordinates of each AGN, measured from our deep Keck images. This
was performed using a custom algorithm (65), which simultaneously models the input sources
at their respective locations, as well as all other nearby sources detected in the SDSS survey
imaging. This modeling fully accounts for the potentially overlapping PSFs of the sources
(a significant issue for WISE given its FWHM^ 6" PSF), and thus effectively de-blends the
photometry. It also utili z es data from multiple visits, thus effectively co-adding all epochs and
improving the effective depth over that in the published catalogs. Table S2 shows photometry
for the five AGN. There we list the SDSS ugriz and WISE 1+1,1+2,1+3,1+4 measurements
determined from our forced photometry procedure, measurements or limits on the Peak 20cm
radio flux F 2 o C m from the FIRST survey, the LRIS R-band photometry, the Lya line-flux, and
the rest-frame Lya equivalent width.
2.1 AGN1: An Obscured Type-2 Quasar
The object AGN1 is embedded in a bright ridge of the nebular Lya emission, as shown in our
Lya image and the 2D spectrum (Fig. 3). Our LRIS spectra have relatively low continuum
signal-to-noise S/N ~ 1 — 2 ratio, consistent with the faint continuum magnitude of AGN1:
V = 23.95 ± 0.07. Nevertheless, we detect strong emission in several lines: HellIA1640,
C ill]A1908, CII]A2328, and a tentative detection of CIVA1549, along with bright Lycc, which
our narrow-band imaging indicated should be strong. The FWHM, line flux, and rest-frame
equivalent widths (EWs) of these lines are listed in Table S3. The presence of emission in
several high-ionization UV lines (i.e. He II and C IV) characteristic of AGN, clearly indicates
that this source is powered by a hard-ionizing spectrum, rather than star-formation. Although
high-ionization lines l ik e C IV, He II, and C III] can be formed in starbursts, stellar winds, the
photospheres of massive stars, and the interstellar medium, the spectra of star-forming galaxies
typically exhibit much smaller rest-frame EWs (< 2A) (66) in these lines than observed in our
spectra. Several independent lines of argument suggest that this object is actually a luminous but
obscured quasar, also referred to as a Type-2 quasar, as we elaborate on below. From the weak
He nA1640 line, we measure z = 2.05476 which is ~ 300km s -1 lower than the estimates from
Lya, and C ill] A1908. We estimate an uncertainty in this redshift of 400km s -1 , which arises
both because of line-centroiding error, and possible systematic uncertainty about the degree to
which a high-ionization line like HellA1640 traces the systemic frame for a Type-2 AGN.
According to unified models of AGN (67, 68), orientation is the primary determining fac¬
tor governing the ultraviolet/optical appearance of an AGN. In this context, vantage points
which are not extincted by the obscuring torus observe a spatially unresolved power-law ul¬
traviolet and optical continuum believed to emerge from the SMBH accretion disk, and broad
high-ionization emission lines (FWHM~ 5000 — 20,000km s -1 ) from photoionized gas in
the so called broad line region at distances ~ 1 pc from the SMBH. However, all vantage
points, observe narrow emission lines from the spatially extended (~ lOOpc — 10 kpc) nar¬
row (FWHM~ 500 — 1500km s -1 ) line region. When the line-of-sight to the central engine is
blocked by a dusty obscuring torus, only the narrow line region is seen. AGN lacking broad
emission lines, but which nevertheless exhibit narrow high-ionization lines are classified as
Type-2 systems, while those with broad-line emission are classified as Type-1. Alternatives to
the AGN unification viewing angle picture argue instead that Type-Is and Type-2s represent
different phases of quasar evolution (69), with all quasars passing through an obscured phase
before outflows expel the obscuring material.
Characteristics of Type-2 quasars at other wavelengths are a hardened X-ray spectrum, re¬
sulting from photoelectric absorption of soft X-rays presumably from the same obscuring torus
extincting the broad-line region, and strong mid-IR (A ~ 10 /mi) emission, which represents
ultraviolet/optical energy emitted from the disk, but absorbed, reprocessed, and re-emitted
isotropically by the torus. In a subset of radio-loud systems, synchrotron radiation, believed
to be emitted isotropically from scales much larger than the torus, is also observed.
Luminous high-redshift radio-loud Type-2 quasars have been studied for decades (see the
review by (70)) but their radio-quiet counterparts have been much harder to identify. Over the
past decade, significant progress has been made in identifying sizable samples of low-redshift
(z < 1) radio-quiet Type-2s from mid-IR (71, 72), X-ray (73), and optical narrow-line emission
(74). However, to date the vast majority of these Type-2s are at z < 1, and there are still
only handfuls of bonafied Type-2 quasars known at z ~ 2 (75) with bolometric luminosities
comparable to the typical SDSS/BOSS Type-1 quasar, L bol > 10 46 ergs -1
The three lines of argument that we use to establish that AGN1 is a Type-2 quasar are 1) the
narrow velocity width of its emission lines 2) the line ratios in its spectrum, and 3) its extremely
red optical-to-mid-IR colors.
Unlike the three other AGN physically associated with the nebula (f/g quasar, AGN2, and
AGN3), which exhibit broad emission lines (FWHM~ 5000 — 10,000km s _1 ), Gaussian fits to
the lines in AGN 1 reveal relatively narrow lines < 1500km s -1 , which are listed in Table S3. At
low-redshifts (z < 0.3), AGN exhibit a bimodal distribution of Ha emission line widths, which
corresponds to the Type-1 and Type2 populations. Based on this distribution, the canonical
dividing line between Type-l/Type-2 quasar classification is FWHM Ha < 1200 km s -1 (76).
However, it is unclear how to extrapolate this classification to the rest-frame UV lines observed
at z ~ 2. First, higher redshift sources are intrinsically much more luminous, have more massive
SMBHs, and hence SMBH-galaxy scaling relations imply that they should be hosted by much
more massive galaxies. Indeed, for the handful of Type-1 quasars at z ~ 2 for which the
narrow-line region [O ill] emission line has been observed in the near-IR, larger line-widths
FWHM= 1300 — 2100km s -1 are indeed observed (77), contrasting with narrower line widths
observed at lower redshift (76). Furthermore, even in the narrow-line region, one generally
expects larger widths for Lya and higher ionization potential lines, because photoionization
modeling generally predicts that emission in these lines is produced at smaller distances from
the nucleus, where velocities are likely to be larger. Indeed, exactly such trends were observed in
the spectrum of a Type-2 quasar at z — 3.288 discovered by (78), who measured FWHM Lya =
1520 ± 30km s _1 , FWHM He ii = 940 ± 140km s -1 , and FWHM cm ] = 1090 ± 140km s -1 ,
which are comparable to our measured values for AGN1. A near-IR spectrum of the Stem et al.
source reveals narrower widths for the rest-frame optical lines that are typically used to classify
low-redshift Type-2s: FWHM^ = 170 ± 130km s -1 and FWHMpm] = 430 ± 30km s -1 .
Note that the unusually hard X-ray spectrum of the Stem et al. source (78) and the implied
high photoelectric absorbing column Ajf = (4.8 ± 2.1) x 10 23 cm 2 leave little doubt that it is
a bonafied Type-2 quasar, notwithstanding the fact that Lya and other high-ionization UV lines
have line widths in excess of FWHM = 1200 km s -1 . We suspect that the canonical rest-frame
optical line-width classification based on z < 0.3 Type-2 quasars (76) needs to be revised for
rest-frame UV lines in more luminous Type-2s at z ~ 2. Based on the narrow lines of AGN1
and their similarity to the Stern et al. source (78), we conclude that it is very likely to be a
Table S3 presents measurements of the line fluxes and measurements or limits on rest-frame
equivalent width for the characteristic AGN lines covered by our spectrum. Based on these
line flux measurements, Fig. S3 shows the line-ratios of AGN 1 in the CIV /He II vs CIv/C III]
plane (left) as well as the C m]/He II vs C III]/C II] plane (right). These are the standard line
ratios conventionally discussed in studies that use photoionization modeling to diagnose the
physical conditions in the narrow-line region (79), although we note that the ratios C III] /He II
vs Cm]/Cll] are not typically plotted against each other. In Fig. S3, we compare the line
ratios of AGN1 to other Type-2 AGN compiled from the literature. Specifically, the circles are
individual measurements of HzRGs (black) from the compilation of (80), and narrow-line X-ray
sources (cyan) and Seyfert-2s (blue) from the compilation of (81). In addition, triangles indicate
measurements from the composite spectra of HzRGs (orange) from (70), the composite Type-2
AGN spectra (purple) from (82), who split their population into two samples above and below
Lya EW of 63A, and a composite spectrum of mid-IR selected Type-2 AGN (green) from (S3).
Finally, the stars represent measurements of these line ratios from composite spectra of Type-1
quasars, based on the analysis of (magenta) (84), (red-magenta) (35), and (blue-magenta) (86).
We caution that determination of robust emission line fluxes for Type-1 quasars is extremely
challenging. Generally the broad emission lines, non-trivial emission line shapes, and the fact
that many of the lines of interest are severely blended, makes the separation of the spectrum
into line and continuum rather ambiguous. To minimize the impact of noise and variations
among quasars, one typically analyzes extremely high signal-to-noise ratio composite spectra,
which also average down quasar-to-quasar variation (84-86). But nevertheless, the resulting
line fluxes and hence line flux ratios are dependent on the method adopted (see (35)) section
4.1 for a detailed discussion). These issues are particularly acute for the lines that we consider:
HellA1640 is heavily blended with the red wing of Civ A1549 and Olll]A1663, AlllA1670,
and a broad Fell complex, Cm]A1909 is severely blended with AlllA 1857 and SiHlA 1892,
and Cll]A2326 is blended with a broad Fell complex. These ambiguities are reflected in the
large variation in the line ratios determined from Type-1 quasar composites in Fig. S3 (stars)
by different studies. Despite these challenges, Fig. S3 clearly illustrates that the line-ratios for
AGN 1 are consistent with the Type-2/Seyfert-2 locus, and differ rather significantly from the
Type-1 line ratios, notwithstanding the somewhat divergent measurements for the latter. We
thus conclude based on the emission line ratios of AGN 1 that it is very likely to be a Type-2
Many studies utilizing WISE data have shown that at bright magnitudes W2 ~ 15 mag,
the WISE Wl — W2 color efficiently selects AGN at z ~ 1 — 3, including obscured objects
(87-89). This results from the rising roughly power law shape of AGN SEDs in the mid-
IR (90,91), and the fact that stellar contamination is low, because the WISE W 1 and W2 bands
cover the Rayleigh-Jeans tail of emission from Galactic stars, resulting in much bluer colors
which cleanly separate from redder AGN. At fainter magnitudes W2 ~ 17, contamination
increases because of optically faint, high-redshift elliptical and Sbc galaxies, which tend to also
have very red Wl — W2 colors. In addition, obscured AGN are also expected to have very
red optical to mid-IR colors, if the mid-IR traces reprocessed dust emission from a heavily
extincted accretion disk emitting in the optical/UV. Indeed, it has been shown that the color-
cut r — W 2 > 6 effectively separates obscured AGN from the their much bluer unobscured
counterparts (89, 92, 93). The WISE photometry in Table S2 indicate that AGN1 is detected in
both W 1 and W2 with W 1 — 16.89 ± 0.08, and W2 = 15.86 ± 0.11, implying a Wl — W2 =
1.03 ± 0.14, thus consistent with the expected red color Wl — W2 > 0.8 of AGN in the WISE
bands. If we adopt our LRIS V-band magnitude as a proxy for r (color corrections are minimal),
then the faint V = 23.95 continuum of AGN1 implies r — W2 = 8.09 ± 0.13, which is also
consistent with the r — W 2 > 6 criterion for obscured AGN. We thus conclude based on mid-IR
and mid-IR to optical colors, that AGN1 is consistent with being an obscured AGN.
2.2 AGN2: A Radio Loud Broad Line AGN
The second quasar (AGN2) is rather faint V = 23.12 (Table SI), but its LRIS spectrum never¬
theless reveals broad Lya, C IVA1549, and C ill] A1909 emission lines, unequivocally indicating
that it is a broad-line Type-1 AGN (see Fig. S2). Centroiding the broad and relatively low S/N
C ill] A1908 line, we estimate ^agn 2 = 2.058 with a 700km s -1 uncertainty. Redshift errors for
rest-frame ultraviolet emission lines are dominated by intrinsic shifts between the emission line
and the true systemic frame, and we adopt the procedure described in (23) to centroid the broad
lines, and use the procedure of (94) to determine the redshift errors.
Motivated, in part, by the previous association of bright Lya nebulae to HzRGs and/or bright
radio sources, we searched the literature and public catalogs for sources in the field surrounding
SDSSJ0841+3921. The Faint Images of the Radio Sky at Twenty cm survey (FIRST) (95)
used the Very Large Array (VLA) to produce a map of the 20 cm (1.4 GHz) sky with a beam
size of 5. "4 and an RMS sensitivity of about 0.15 mJy beam -1 . The survey covers the same
10,000 deg 2 sky region covered by the SDSS imaging, has a typical detection threshold of
1 mJy, and an astrometric accuracy of 0.05". The FIRST catalog reports only a single source
within 60" of SDSSJ0841+3921 at RA=08:41:58.6 DEC=+39:21:14.7, approximately 6" South
of the f/g quasar and coincident with AGN2. The FIRST peak and integrated fluxes at 20 cm are
Fpeak = 12.87 mJy/beam and F int =14.81 mJy respectively, giving a morphological parameter
0 = log 10 (F int /Fp eak ) = 0.06. This indicates a compact, unresolved source at 1.4 GHz (96).
Therefore, there is no evidence for a bright, extended radio source (e.g. a radio jet) associated
with the protocluster system and giant Lya nebula.
Follow-up VLA observations of this source were performed by (97), having erroneously as¬
sociated it with the f/g quasar. The spectral slope measured from 8.4 GHz and 4.9 GHz observa¬
tions is relatively steep, af g = —1.21, which is generally interpreted as evidence that the source
has a large viewing angle (i.e. edge-on) (98). This may indicate that FIRSTJ084158.6+392114.7
(AGN2) is oriented away from us, although we re-emphasize that it exhibits a Type-1 spectrum.
These much deeper radio images of the source show no evidence for extended structure or jets,
providing further evidence that it is not a HzRG (97).
2.3 AGN3: A Faint Broad Line AGN
Similar to AGN2, the source AGN3 is also faint V = 23.09 but exhibits broad emission-lines of
CIVA 1549, C III] A 1908, and Mg ilA 2798 (see Fig. S2). Unfortunately Lyct was not covered by
our DEIMOS spectrum, but the source has a large rest-frame equivalent width (VF Lya = 35A)
indicative of strong emission. Based on the broad lines, we also classify this source as a broad¬
line (Type-1) AGN. We determine a redshift of zagn 3 = 2.050 from the C ill] A1908 line with a
700km s' 1 uncertainty.
3 Probability of Finding a Quadruple Quasar
The probability dP of detecting three quasars around a single known quasar (here the f/g quasar)
can be written as
dP = n% so dV 2 dV 3 dV 4 [1
T£(?T 2 ) + --- (6 permutations)
+ C(' r i 2 ,r’ 23 ,r' 3 i) H- (4 permutations) (2)
+ (((r* 12 )^ 34 ) + ••• (3 permutations)
+ v( r lU 2U3U4)]
where tiqso is the number density of quasars, dV, is the infinitesimal volume element centered
on the location r* of the 2 th quasar, is the distance between the /th and jth quasar, and £,
and r) are the two-point, three-point, and four-point correlation functions, respectively (99).
The total probability P of finding a quadruple quasar system within a maximum radius r max
is the integral f dP over the three volume elements dV r , for all possible configurations of the f t .
Our goal in what follows is to obtain an order of magnitude estimate for this total probability.
First note that on the small scales of interest to us (here r < 200 kpc), many of the terms in
eqn. () can be neglected. For the higher order correlation functions ( and ?/, a scaling of the
form ( oc £ 2 and ?/ oc £ 3 is often assumed. This arises from the fact that for Gaussian initial
conditions and a scale-free initial power spectrum, this scaling is obeyed to second order in
Eulerian perturbation theory (100). However, this scaling is not valid for the highly non-linear
small-scales of interest to us here. On these small scales, we follow the halo model approach
(101), and write the higher order functions as a sum of terms representing contributions from
the multiple possible halos. For example the four-point function can be written as the sum
V = + n 2h + rf h + 27 4h , (3)
resulting from contributions from one to four halos. As we are concerned with scales com¬
parable to the virial radius of the dark matter halo hosting, i.e. the f/g quasar r < 200 kpc,
the one halo term, which quantifies the four-point correlations when all four points lie in the
same dark matter halo, will dominate (101). It is easier to calculate this one-halo term of the
four-point function ?/ lh in Fourier space, where one instead works with the one halo term of the
trispectrum T lh , and it can be shown that T lh scales as the Fourier transform of the dark matter
halo density profile to the fourth power T lh oc p^ M . This scaling holds likewise for r/. Thus for
quasar 4-point correlations on small-scales, the halo model indicates that we expect 77 oc Pq SO ,
where Pqso is the density profile of quasars in the dark matter halos which host them.
By a similar line of argument, the small-scale two- and three-point functions scale as £ oc
Pq SO and (' oc Pq SO , respectively. Or alternatively, (' oc £ 3 / 2 and r/ oc £ 2 . On the proper scales
of interest to us here r < 200 kpc (r com = 600kpc, comoving), £ 1 and thus the terms
dominating eqn. () are the three permutations of terms like £(r 12 )£(r 34 ), and //, which are all
of comparable order £ 2 .
Thus in order to compute the total probability of finding three quasars around a known
quasar, we must integrate the four dominant terms in eqn. Q, over all possible configurations
of the Tj, which can be realized within a spherical volume, whose radius is set by the largest
separation in the quadruple quasar r max . These integrals all have a common form and order-
of-magnitude. For example the first three terms involving pairs of correlation functions can be
£(r)47rr 2 dr
C(r M )dV :i dV 4 ~ n\
£(r)47rr 2 c/r ) V, (4)
where V = 47r/3rf nax . Thus the final expression for an order-of-magnitude estimate of the
( pr max \ 2
J 1 £(r )4:7Tr 2 drJ V , (5)
where the factor of four accounts for the four similar dominant terms of order ~ £ 2 in eqn. (|2j).
To calculate n qso, we use the quasar luminosity function of (102) assuming an apparent
magnitude limit V < 24. This limit is highly conservative, as even AGN2 and AGN3 have
V = 23.1 mag. Although the Type-2 AGN1 has V = 24, its nucleus is obscured from our
perspective, and this magnitude represents host-galaxy emission or a small-fraction ~ 5% of
scattered nuclear emission. The corresponding nuclear V -band emission of this AGN is likely
to be much brighter than V = 24 given that its Lyo: flux, coming from the narrow-line region, is
so strong, i.e. larger than that of AGN2. To account for the fact that our luminosity function only
includes unobscured AGN, we divide it by a factor (1 — f Qbs ), where / obs is the obscured fraction
of AGN. We adopt an obscured fraction of f Q bs = 0.5 consistent with recent determinations
(103,104). Altogether, we estimate that the comoving number density of faint AGN is 71qso =
3.8 x 10 -5 cMpc~ 3 .
We assume that all four AGN reside in a physical structure with physical size r max . This is
a good assumption, given that: the projected separations of the AGN are small, they reside in a
larger-scale overdensity of Lycc emitters, the f/g quasar, AGN1, and AGN2 are all embedded in
the giant Lycc nebula, and the morphology of the nebular emission around the AGN suggests a
physical association. Moreover, although we argue here that the probability of finding a quadru¬
ple AGN is very small, it would be much smaller if the AGN had larger line-of-sight separations,
in which case one cannot invoke as large of an enhancement due to clustering. Taking the f/g
quasar location as the origin, the largest angular separation measured in the quadruple quasar is
that of AGN3 with 6 = 18", corresponding to to a transverse separation r± = 150 kpc. We thus
choose r max = 250 kpc, such that the the spherical comoving volume we consider corresponds
to V = 1.9 cMpc 3 . The probability of finding the three AGN around the f/g quasar at random,
i.e. in the absence of clustering, is extremely small P ran = (nQso ^) 3 = 3.4 x 10 -13 .
The small-scale two-point correlation function of quasars was first computed by ( 6 ), and
later extended to even smaller scales ~ 5 kpc by the gravitational lensing search of (30), who
found that a power law form £ = (r/r 0 ) -7 with 7 = —2 and r 0 = 5.4 ± 0.3 /r _1 cMpc provides
a good fit to the data over a large range of scales. Plugging these numbers into eqn. 0, we
finally arrive at a probability of P = 6.9 x 10 -8 , justifying our order of magnitude estimate of
P ~ lx 10~ 7 .
4 Giant Nebulae-LAE Clustering Analysis
4.1 Constructing the LAE Catalog
We created an LAE catalog from our images of SDSSJ0841+3921 using SExtractor (105) in
dual-mode, using the NB image as the reference image. In order to minimize spurious detection
we varied the parameters DETECT JYIINAREA (minimum detection area) and
DETECT.TIIRESII (relative detection threshold) creating a large number of LAE candidate
catalogs. We verify the number of spurious detections for each catalog running SExtractor on
the “negative” NB image obtained by multiplying the NB image by —1. The ratio between the
number of sources detected in the NB image and the sum of sources detected in both the NB
and “negative” image define a “reliability” criterion. We selected DETECT_MINAREA equal
to 8 pixels and DETECT.THRESH equal to 1.8cr corresponding to a ’’reliability” of 94%. We
used the same parameters to obtain a catalog of sources for the E-band image. In order to select
LAE candidates, we estimated the NB — V color from SExtractor isophotal magnitudes and
applied a color-cut corresponding to a Lya rest-frame equivalent width of W\, yoi > 20A which
is the standard threshold used in the literature. In particular, we assumed a flat continuum slope
in frequency to convert the V-band magnitude to a continuum flux at the Lya wavelength. For
the objects in the catalog without broad-band detection (i.e., with S/N< 3), we determine a
lower limit on the EW following (61). Finally, we inspected each object by eye and removed
four spurious sources that were lying in proximity to a bright star.
The final clean catalog consisted of 61 LAE candidates above a Lya flux of limit of 5.4x 1CT 18
erg s _1 cm -2 , over our 6 ' x 7.8' imaging field-of-view, and the redshift range z = 2.030 — 2.057
spanned by the narrow-band filter. We estimated that our catalog is 50% (90%) complete above
a flux limit of P 50 = 6.7 x 10 -18 erg s _1 cm -2 (P 90 = 7.4 x 10 -18 erg s _1 cm -2 ). This catalog
of LAEs is presented in Table S4. In the following clustering analysis, we consider only those
sources with rest-frame Wj Aya > 20A and a luminosity log 10 L^a > 42.1 which at z = 2.04
corresponds to a flux limit of 4.0 x 10 — 1 ' erg s _1 cm -2 . Given that this flux level is a factor
of five higher than the 90% completeness flux limit of our catalog, we are certain that we are
complete to such sources. The f/g quasar lies at a distance of 2.08' from the edge of our imag¬
ing field, and for simplicity, our clustering analysis considers only the 10 sources within this
radius which corresponds to a comoving impact parameter of 2.2 h~ l cMpc. Of these, 3/10 are
the spectroscopically confirmed AGN1-3. The sources included in our clustering analysis are
denoted by an asterisk in Table S4.
4.2 Clustering Analysis of HzRGs and LABs
We now describe how we used a compilation of LAEs around HzRGs and LABs to estimate
the giant nebulae-LAE cross correlation function. The vast majority of z ~ 2 — 3 protoclusters
in the literature have been identified and studied via overdensities of LAEs over comparable
fields of view 6' x 7.8' as our LRIS observations. For the HzRGs, we used a catalog of LAE
positions around HzRGs from the survey of (13). Specifically, we focus on six HzRGs from
the Venemans et al. study in the redshift range z = 2.06 — 3.13, which are approximately
co-eval with SDSSJ0841+3921 at z — 2.0412. In addition to the six HzRGs from Venemans
et al. (13), we also include two LABs at z ~ 2 — 3 from the literature whose environments
have been surveyed for LAEs. Nestor et al. (106) conducted a narrow band imaging survey
of the well known SSA22 protocluster field at z — 3.10. Their Keck/LRIS observations were
centered on LABI, whereas LAB2, the other known LAB in this field, resides at the edge of
their imaging field. For this reason we only consider the overdensity around LAB 1. Prescott et
al. (19) conducted an intermediate-band imaging survey for LAEs around an LAB at 2 = 2.66
(LABd05) discovered by Dey et al. (16), finding a significant overdensity, and argued that the
LAB resides in a protocluster. The eight objects used in our clustering analysis are listed in
Given that the Venemans et al. sample (13) is the largest compilation of LAEs around
HzRGs, and that, to our knowlege, the two LABs we consider are the only ones whose envi¬
ronments have been characterized using LAEs to a depth comparable to our observations, this
sample of eight HzRGs/LABs is the only dataset currently avaailable for conducting our clus¬
tering analysis and we believe that they comprise a fairly representative sample. Venemans et
al. applied standard HzRG selection criteria to define their sample, and they published results
for all HzRGs, not just those that hosted dramatic overdensities. Although LABI and LABd05
are larger and brighter than more typical LABs, they are comparable in size and luminosity to
the HzRGs and to SDSSJ0841+3921, making for a fair comparison. One cause for concern is
that the distribution of LAEs around the LABs were published specifically because these two
systems resided in overdensities, and thus there may be a publication bias against LABs resid¬
ing in lower density environments. Given that LABS only comprise a fourth of our clustering
sample, we do not believe that this significantly biases our results, but repeating our clustering
analysis for a larger and more uniformly selected sample of LABs would clearly be desirable.
Note that the HzRG MRC 1138—2262 included in our analysis is the famous Spider Web
Galaxy protocluster, which has been the subject of extensive multi-wavelength follow-up, and
is one of the most dramatic protocluster systems known. Venemans et al. (13) argued that all of
the HzRGs in Table S5 reside in overdensities of LAEs, with the exception of MRC 2048—272
and TN J20093040, for which the abundance of LAEs was found to be consistent with the field
value, albeit with large error bars. We nevertheless include these two sources in our clustering
analysis for several reasons. Our objective for the clustering analysis is to shed light on the con¬
nection between active SMBHs, large-scale Lya nebulae, and protocluster environments. Given
that it has been argued that the HzRGs as a population trace extremely overdense protocluster
environments, and given that MRC 2048272 and TN J20093040 are both powerful HzRGs with
large Lya nebulae, excluding them from the analysis would be arbitrary, and would bias our
correlation function high. Furthermore, the background number density of LAEs used by Vene¬
mans et al. (13) is rather outdated, the overdensity estimates have large error bars for individual
systems, and the luminosity limits that they adopted for the LAEs in each field considered
were heterogeneous. These factors complicate the interpretation, particularly if the clustering
of LAEs is luminosity dependent. Hence our goal is to conduct a homogeneous analysis of
the clustering of LAEs around the HzRGs/LABs which are above a uniform Lya luminosity
log 10 Lhya > 42.1, without cherry picking specific objects. We adopt this luminosity limit for
two reasons. First, the NB imaging observations for all eight objects listed in Table S5 are com¬
plete for identifying LAEs above this luminosity, mitigating uncertainties due to incomplete¬
ness. Second, this limit corresponds to the faintest luminosity for which the the background
number density of LAEs is well characterized (32); working fainter would thus require a signif¬
icant extrapolation of the luminosity function, which would add systematic uncertainties to our
results. For all of the objects analyzed (see Table S5), the LAE catalogs available are complete
above the Lya equivalent width limit of WL ya > 20A (the same value used to define our LAE
catalog around SDSSJ0841+3921 ), with the exception of the Prescott et al. observations of
the LAB LABd05. These observations used an intermediate band filter, and thus adopt a rest-
frame W\, ya > 40A. We consistently account for this difference below when we compute the
background number density of LAEs.
Given the small size of this HzRG/LAB dataset, we use an unbinned maximum likeli¬
hood estimator to determine the parameters r 0 and 7 , following the procedure outlined in
(107,108). Specifically, we assume that the cross-correlation function between the giant nebu¬
lae (HzRGs/LABs) and LAEs obeys a power law form £(r) = (r/r 0 ) -7 . We use this correlation
function to calculate the expected number of LAEs within a comoving cylindrical volume with
transverse separation from R to R + dR, and half-height A Z, which is set by the width of
the given narrow band filter A Z = FW 11X1/ (2a 11 (z)), where FWHM = cAA/A is the full
width half maximum of the filter in velocity units, a = 1/(1 + z) is the scale factor, H(z ) is
the Hubble expansion rate, and division by aH(z ) converts velocities to comoving units. We
imagine dividing the transverse separation R into a set of infinitesimal bins of width dR, such
that each bin can contain only one or zero LAEs. Under the assumption that the clustering of
LAEs around the giant nebulae is a Poisson process, we can write the likelihood function as
( 6 )
where /i = 47rf?AZn L AE[l + h(R)]dR is the probability of finding a pair in the interval d,R, the
index i runs over the dR bins containing the N L ae HzRG/LAB-LAE pairs in the sample, and
the index j runs over all the bins for which there are no pairs. Here tilae(> Lhya) is the number
density of LAEs brighter than the survey limit L^ ya , and h(Rj is the cross-correlation function
averaged over the filter width
Taking the natural logarithm of the likelihood, we have
( 8 )
where we have dropped all additive terms which are independent of our model parameters
(r 0 ,7), and [ /f. tr , in ,7T. m;1X J is the range of comoving scales over which we search for HzRG/LAB-
LAE pairs. In the following clustering analysis, we adopt R m in = 50 h~ 1 ckpc to minimize
confusion with LAEs embedded in the respective nebulae, and R max = 2.5 h~ x cMpc, which
is well matched to the maximum impact parameter probed R = 2.2 h~ l cMpc around the f/g
quasar. Given that the HzRGs and LABs in Table S5 all have different filter widths, as well as
different VE Lya limits and redshifts (which changes the background density tilaeX the likelihood
in eqn. <) is specific to a single source. However, the full likelihood of the entire dataset given
the model parameters, can be determined by simply summing log-likelihoods (multiplying the
likelihoods) over all of the HzRGs/LABs, each of which has the form of eqn. (),
Our clustering analysis relies on a determination of the field LAE luminosity function tilae-
We use the Ciardullo et al. (32) Schechter function fits to the luminosity function from their own
LAE survey at z = 3.1, and the Ciardullo et al. fits to the LAE survey data of (109) at z = 2.1.
Both of these surveys are highly complete down to our limiting luminosity of log 10 L hya > 42.1
for LAEs with equivalent width W^a > 20A. Ciardullo et al. also estimated the rest-frame
equivalent width distribution, and found that LAEs follow an exponential distribution with rest-
frame scale lengths of W 0 = 50A and W 0 = 64A, at z = 2.1 and z — 3.1 respectively.
Hence, for LABd05 which employed an W limit > 40A, we take the luminosity function of
LAEs to be n LAE (> W hmit , > L hya ) = exp [-(IE limit - 20A)/IE 0 ]^lae(> 20A, > L LyQ ),
where tilae(> 20A, > L^ ya ) are the Ciardullo et al. Schechter-function fits. To determine
the luminosity function at intermediate redshifts between z = 2.1 and z = 3.1, we linearly
In Fig. S4 we show confidence regions in the r 0 — 7 plane for the likelihood in eqn. ([ 8 ]). The
maximum likelihood value is r 0 = 29.3 h~ l cMpc and 7 = 1.5, however a strong degeneracy
exists between r 0 and 7 , such that neither is individually well-determined. The low signal-to-
noise ratio of the clustering measurement is partly to blame for this degeneracy, however it is
also a consequence of the functional form adopted for the correlation function £(r) = (r/ro) -7 ,
for which the amplitude and slope of the correlation function are not independent parameters.
Note that errors on the parameters r 0 and 7 given by our maximum likelihood estimator
in eqn. ([ 8 ]) and illustrated in Fig. S4 underestimate the true errors. This is because our likeli¬
hood explicitly assumes that the positions of the LAEs are completely independent draws from
a Poisson process. This approach ignores correlations between the LAE positions (which are
generically present for hierarchical clustering) and fluctuations due to cosmic variance. We
opted for this unbinned estimator, because the traditional method of computing the correlation
function in impact parameter bins, and fitting a power-law will, for such small samples, yield
results which depend sensitively on the choice of binning, unless the covariance of the bins is
included in the fit. The correct approach would be to fit the binned correlation function using the
correct covariance matrix, whose off-diagonal elements would reflect the correlations between
LAEs, and whose amplitude would include fluctuations due to sample variance. Given our lim¬
ited sample size of eight HzRGs/LABs, the total number of LAEs in our clustering analysis is
only 75, and there is no hope of computing the covariance from the data itself. In principle one
could adopt a model for the covariance based on IV-body simulations, but this would require
a full model for how LAEs populate dark matter halos, and would be sensitive to the assumed
masses of the protocluster halos. Such a detailed analysis is beyond the scope of the present
work. As a compromise, we use the unbinned maximum likelihood estimator in eqn. © to
determine our best fit model parameters, and we obtain estimates of the errors on parameters
via a bootstrap technique. Specifically, we generate 1000 realizations of mock datasets, where
a random sampling of the eight HzRGs/LABs in Table S5 are selected with replacement. Given
that we are measuring a cross-correlation of LAEs around eight distinct objects, only the LAE
positions arising from the same source are correlated, whereas the others are completely in¬
dependent because of the large respective distances between the HzRGs/LABs. As such, for
each of these 1000 realizations, we also bootstrap resample the positions of the LAEs around
each object (with replacement), which encapsulates fluctuations due to correlations in the LAE
positions and Poisson counting fluctuations. Lor each bootstrap realizationn of our dataset, we
compute the parameters r 0 and 7 via our maximum likelihood procedure, and determine errors
on these parameters from the resulting distributions.
Based on this procedure, we define the lex confidence region from the 16th and 84th per¬
centiles of the bootstrapped parameter distributions. We find a range 0.36 < 7 < 1.98 for 7 ,
and 15.8h _1 cMpc < r 0 < 9404 h -1 cMpc. The correlation length r 0 is poorly constrained,
because for flat correlation function slopes 7 < 1 , the clustering projected along the line-of-
sight is very insensitive to r 0 , which again simply reflects the fact that the amplitude and slope
of the correlation function are not independent parameters. This is the same behavior seen in
the contours in Lig. S4. Lixing the slope 7 to its maximum likelihood value of 7 = 1.5, we
determine a la confidence region for the cross-correlation length of r 0 = 29.3 ± 4.9 hr 1 cMpc.
Although we do not fit the binned correlation function to determine model parameters, it is
instructive to compute a binned correlation function to visualize the data, our model fits, and
the errors on these fits. As such, we define a dimensionless correlation function as
X(Rii Ri+1 )
- 1 ,
where (PG)ij is the number of real HzRG/LAB-LAE pairs in the impact parameter bin \R t , R i+ 1 ]
around the j th object. The quantity (PR)ij = U],ae(Rj, > Ibiimitj)^' i- s the expected random
number of pairs, which is simply the product of the background number density of LAEs, mlae,
and the cylindrical volume, 14, , defined by the bin [ R,, R i+ x ] and the half-height AZj of the NB
filter used to survey the jth object. We compute errors on the binned correlation function via
an analogous bootstrap method, whereby x(R-n Ri+i) is recomputed from 1000 bootstrap re¬
samplings of the sample of eight HzRGs/LABs while simultaneously bootstrap resampling the
LAE separations from each object.
In Fig. S5, we show the binned correlation function calculated via this procedure, with la
error bars, determined by taking the 16th and 84th percentiles of the distribution of x(Ri, Ri+i)
resulting from the bootstrap samples. For any value of the model parameters r 0 and 7 , we can
compute the predicted value of y R i+ \) by evaluating the integrals
rAZj rRi+ 1
(PG)ij = n LAE ( Zj , > Wi imitJ )I4 / / [1 + aVR 2 + Z 2 )}2 -kR dRdZ, (10)
and then evaluating the sum in eqn. fii. The uncertainty on the predicted clustering level,
arising from the uncertainty on our model parameters, can thus be determined by evaluating the
16th to 84th percentile confidence region of the predicted x(Ri, R;+\ ), using the distribution of
r 0 and 7 from the bootstrap resampling of our model fits. The resulting lcr uncertainty on the
predicted correlation function is shown as the gray shaded region in Fig. S5, with the maximum
likelihood value r 0 = 29.3 h _1 cMpc and 7 = 1.5 shown by the red curve. We see that our
maximum likelihood model fit and bootstrap errors are consistent with the binned correlation
function estimates and their respective errors.
Given the results of the clustering analysis described in this sub-section, we generate Fig. 2
of the main text as follows. Using our catalog of LAEs around the f/g quasar with log 10 L Lya >
42.1 (Table S4), and the LAE luminosity function prediction for the background number density
^lae, we can compute the cumulative overdensity profile 5(< R) = 1 + £(-R min , R) = ( PG)/
(PR ). The overdensity profile predicted by our maximum likelihood determination of r 0 and
7 for the giant nebulae-LAE correlation function are determined by evaluating eqn. (10). This
is shown as the red curve in Fig. 2. Similarly we can determine the lcr error on the predicted
d(< R) by evaluating it for all of our bootstrap samples of r 0 and 7 and taking the 16th-84th
percentile region of the resulting d(< R) distribution. This region is shown as the shaded grey
area in Fig. 2. For comparison we show the expected overdensity (projected over our NB filter)
of Lyman break galaxies (LBGs) around radio-quiet quasars, based on the recent measurements
of quasar-LBG clustering at z ~ 2.7 by Trainor et al. (JO), who found a cross-correlation
length of r 0 = 7.3 ± 1.3 /r _1 cMpc, assuming a fixed slope of 7 = 1.5, similar to the best-fit
slope deduced for our giant nebulae-LAE correlation. The blue dotted lines indicate the la
error quoted by Trainor et al. (10) on the cross-correlation length with the slope held fixed.
Comparing to the Trainor et al. (10) measurements may seem questionable, given that: a) this
is the cross-correlation with LBGs, whereas the clustering analysis in Fig. 2 is concerned with
LAEs; b) the quasars studied by Trainor et al. were hyper-luminous. However, both of these
differences are likely to produce small effects on the cross-correlation strength because a) the
clustering of LAEs and LBGs are comparable at z ~ 2.7, and, assuming hierarchical clustering,
the cross-correlation strength is only sensitive to the square-root of any differences in clustering
strength; and b) the Trainor et al. (10) results are consistent, within the errors, with the quasar-
LBG cross-correlation results of (110), who studied much fainter quasars. Indeed, both (110)
and Trainor et al. (10) argue that the clustering of LBGs around quasars is independent of
luminosity, assuaging concerns about comparing to hyper-luminous quasar clustering.
5 Near-IR Observations: The Systemic Redshift of the f/g
Quasar redshifts determined from the rest-frame ultraviolet emission lines (redshifted into the
optical at z ~ 2 ) can differ by up to one thousand kilometers per second from the systemic
frame, because of the complex motions of material in the broad line regions of quasars (94,111).
To achieve higher precision, one may analyze spectra of the narrow [O in] emission lines emerg¬
ing from the narrow-line region, but at z ~ 2, this requires near-IR spectroscopy. On UT 2007
April 24, we obtained a near-IR spectrum of the f/g quasar (and also the b/g quasar) with the
Near InfraRed Imager and Spectrometer (NIRI) on the Gemini North telescope. These data
were reduced and calibrated with the LowRedux package, which performs wavelength calibra¬
tion using the atmospheric sky lines recorded in each spectrum. With the chosen configuration,
the spectra span 1.40 — 1.94//m in the H-band.
Fig. S 6 presents the NIRI spectra, fluxed with telluric standard stars taken that night. Each
quasar shows the detection of strong [O ill] emission and modest H (3 emission. The [O ill] A5007
lines were centroided using a flux-weighted line-centering algorithm with pixel values weighted
by a Gaussian kernel, and this algorithm was iterated until the line center converged to within
a specified tolerance (see (25) for more details). A small offset of 27km s -1 was also added
to the redshifts determined in this manner to account for average blueshift of the [Olll]A5007
from systemic (111). We measured wavelength centroids of A 0 b s = 1.5231 ± 0.0001/rm and
1.6095 ± 0.00006// 111 for the f/g and b/g quasars, respectively. This gives systemic redshifts
% = 2.0412 for the f/g quasar and 2.2138 for the b/g quasar. We estimate an uncertainty in
this redshift of 50 km s' 1 based on error in our wavelength calibration, line centering, and the
estimated uncertainty when one uses [O ill] to assess the systemic redshift of a quasar ( 111 ).
In the Main text, we described the complex kinematic motions of the quadruple AGN system
and nebula as measured with slit spectroscopy (optical and near-IR; Figs. 3,S6) and absorption¬
line spectroscopy of a b/g quasar (Fig. 4). The AGN exhibit line-of-sight velocity differences
of at least 8v = 1300km s _1 , requiring extreme gravitational motions. Note that these large
velocities cannot be explained by Hubble expansion. The miniscule probability of finding a
quadruple quasar system in the absence of clustering P ~ 10 -13 (see Supp.[|]), and the physical
association between the AGN and giant nebula, demand that the four AGN reside in a real
collapsed structure, and thus the relative motions of the AGN must result from gravitational
motions in a collapsed structure.
The velocities of the nebular Lya emission are also extreme, exhibiting line-of-sight veloc¬
ities relative to the systemic redshift of f/g quasar ranging from 5v = —800 to +2500 km s -1 .
These also trace the gravitational potential of the protocluster but may be influenced by other
physical processes (e.g. kinematic or radiative feedback from the AGN). Finally, the absorp¬
tion line kinematics of cool metal-enriched gas in the protocluster halo, measured from the b/g
quasar spectrum, show strong absorption at ~ +650 km s -1 with a significant tail to velocities
as large as ~ 1000 km s -1 .
We further emphasize two aspects here. First, although the nebular emission exhibits large,
relative velocity offsets, the internal motions are relatively quiescent and generally unresolved
at the Keck/LRIS spectral resolution (FWHM « 160km s -1 ). There is no evidence for the
“double-peaked” emission characteristic of resonantly trapped Lya. This implies that the gas
within the nebula has modest opacity to Lya photons and that the observed radiation is not
predominantly scattered photons. These data also indicate that the gas motions are not strictly
random, but are rather organized into coherent flows, i.e. the velocity offsets observed greatly
exceed the measured velocity dispersions. Similar coherent motions have been previously ob¬
served in giant Lya nebulae in HzRGs ( 112-115 ) and LABs (116).
Second, consider the gas kinematics along the single Slit 2 (Figs. 3and SI). This slit covers
the nebula from near the f/g quasar, where one observes a negligible velocity offset (8v ~
0km s -1 ), to the b/g quasar at « 200 kpc separation. As one travels along the slit, the velocity
offset increases to ~ +500km s -1 as one approaches the b/g quasar, where the absorption-line
spectrum shows a velocity offset of 8v ~ +650km s” 1 (Fig. 4). Altogether, the gas along
Slit 2 exhibits a velocity gradient of at least 500km s -1 across 200kpc (projected). These data
establish the presence of large-scale, high-velocity flows of cool gas within this protocluster
7 Estimates for the Abundance of Giant Lya Nebulae Around
The protocluster SDSSJ0841+3921 was discovered because of its large-scale Lya emission
properties, and was selected from a well-defined survey for extended Lya emission around 29
quasars (22). We can use the statistics of this survey to estimate the prevalence of giant Lya
nebulae around quasars. Insofar as there appears to be a physical connection between giant Lya
nebulae and protoclusters, this also provides insight into the prevalence of protoclusters around
quasars. We also present an estimate for the frequency of occurrence of giant Lya nebulae from
an independent narrow-band imaging survey of quasars conducted by our group, and also dis¬
cuss the results of similar narrow-band imaging surveys performed by others. Although much
past work has been dedicated to characterizing Lya nebulae around radio-loud quasar sam¬
ples (11 7), it is only recently that sensitive observations have been undertaken to characterize
the abundance of nebulae around predominantly radio-quiet quasars.
Our slit-spectroscopic observations allow us to estimate the effective covering factor of ex¬
tended Lya emission around the typical quasars (L bo i ~ 10 46 ergs' ') we surveyed. Follow¬
ing the discussion in § 4.3 of (22; see Fig. 4 of that paper), large ~ 70 kpc scale nebulae
can be detected down to 2 x SB !rT , where SB lfJ are the la surface brightness limits quoted
in Table 1 of (22). In that sample a nebula of comparable brightness to SDSSJ 0841+3921
SB Lyo , ~ 10 - 17 ergs -1 cm -2 arcsec -2 could have been detected at > 2 a significance in 23/29
cases (see Table 1 of (22)). Based on this, we can crudely determine the covering factor of
such filamentary emission as being the ratio of the area of the filament covered by our origi¬
nal slit observation of SDSSJ0841+3921 (Slit2 in Figs. 3and SI), to the total area surveyed at
sensitivity sufficient to detect a SBL ya > 10 -1 ' ergs -1 cm -2 arcsec -2 nebula.
We consider an annular region around a z = 2 quasar with an inner radius of R ±An \ n =
50 kpc (5.8") and outer radius of R L , max = 150 kpc (17.3"). The inner radius is chosen to ex¬
clude the small-scale fuzz at R± < 50 kpc, which is detected in a large fraction of quasars (22),
from the larger scale filamentary emission. The outer radius corresponds to the maximum extent
of the (high significance) emission detected in our discovery spectrum of SDSSJ 0841+3921
(see Slit2 in Fig. 3panel c, and the middle column of Fig. SI), and is about the virial radius of
the dark matter halos M vir = 10 12,5 M 0 at z = 2 which host quasars (11). We also observe
a drop-off in the covering factor of optically thick absorbers at R± ~ 200 kpc, comparable
to this R gmax (27). We assume the filamentary emission in the discovery spectrum extends
from +50 kpc to +150 kpc at positive slit position, and —50 to —80 kpc at negative slit po¬
sitions, corresponding to an area of 15 sq. arcsec given our 1" slit (1130 kpc 2 ). The total area
is equal to the Nq SO = 23 quasars observed to sufficient depth, times the 23 sq. arcsec area
(1740 kpc 2 ) for which our slit intersects the annular region, implying a total area of 530 sq.
arcsec (39922 kpc 2 ) surveyed. For reference, the total annular area around a single quasar is
834 sq. arcsec (62832 kpc 2 ), so the effective area of our survey is 64% of the area around a
single quasar. Hence we deduce a covering factor of filamentary emission of f c ,m = 0.028.
Thus taken in aggregate, we find that only ~ 3% of the area between R± = 50 — 150 kpc
around quasars is covered by Lya emission with SB LyQ . > 1CT 17 erg s -1 cm -2 arcsec -2 , which
can be compared to the much larger covering factor of 31% above this level, measured for
emission over the same area from our NB image of SDSSJ0841+3921 (see Fig. la). The ag¬
gregate and individual covering factors can be commensurated if the majority of quasars do
not exhibit large-scale Lya nebulae above SB Lya > 10 - 17 ergs -1 cm -2 arcsec -2 , but a small
fraction ~ 10% have bright nebulae covering a large area, as in SDSSJ0841+3921, and the
product of this ~ 10% frequency and ~ 31% covering factor set the aggregate ~ 3% effective
covering factor for the entire population. This argument is consistent with our intuition that,
given the random orientation of the slit, we could probably only have identified the asymmetric
extended emission around SDSSJ0841+3921 for ~ 50% of the possible slit orientations, and
given that 23 objects were surveyed to the required depth, this again suggests that ~ 10% of
quasars exhibit giant nebulae.
An independent estimate for the frequency of giant Lya nebulae around quasars can also be
determined from a narrow band imaging survey conducted by our group. As part of a continuing
program to search for fluorescent Lya emission powered by luminous z ~ 2 quasars, to date we
have obtained deep narrow band imaging of a sample of 26 objects. These survey observations
are analogous to the narrow-band imaging of SDSSJ 0841+3921 discussed in Supp. 1.2, in that
they use custom designed narrow-band filters, however the targets constitute an independent
sample of quasars, which are not part of the QPQ survey of (22). These observations have been
conducted on multiple telescopes, which we briefly summarize. For the present purposes we
define a giant Lya nebula to be extended emission on a scale > 50 kpc from the quasar with
average SB Lya > 10 -17 ergs -1 cm -2 arcsec -2 .
The quasar HE0109—3518 at z = 2.41 was observed with the Very Large Telescope/ -FOcal
Reducer Spectrograph (VLT-FORS) using a custom designed narrow band filter (A cent er =
4145A, FWHM a = 40A) for a total exposure time of about 20 hours, achieving a sensitiv¬
ity limit (la) of SBL ya = 8 x 10 -19 ergs -1 cm -2 arcsec -2 (61). No extended Lya nebula
was discovered. We have imaged a total of eight quasars with the Keck LRIS imaging spec¬
trometer, of which three were observed with a filter centered on Lya redshifted to z — 2.279
(Acenter = 3986A, FWHM A = 30A), and five with one centered at z = 2.190 (A cen ter = 3878A,
FWHM a = 30A). The exposure times for the LRIS observations vary from 2 to 10 hours, re¬
sulting in sensitivity limits (la) ranging from SBL ya = 1.2 x 10 -18 ergs -1 cm -2 arcsec -2 to
5 x 10 -19 . These VLT/LORS and Keck/LRIS observations specifically targeted hyper-luminous
quasars V < 17, with the goal of understanding fluorescent Lyct emission from dark-galaxies
(61) and the quasar CGM (41), and led to the discovery of a giant Lya nebula around the
quasar UM287 at z — 2.28. Nebulae were not detected around any of the other quasars
surveyed. Linally, using a custom filter centered on Lya at z — 2.253 (A cen ter = 3955A,
FWHM a = 30A), we have observed a total of 17 quasars using the Gemini Multi Object
Spectrograph (GMOS) on the telescope Gemini South. For the GMOS observations, a subset
of 3 quasars were observed with longer integrations, typically 5 hours and achieving a depth
of (la) SB Lyo , = 1.5 x 10 - 18 ergs -1 cm - 2 arcsec -2 ; whereas the other 14 quasars were ob¬
served in a fast survey mode with typical exposure times of about 2 hours and depth SB Lya =
3 — 4 x 10 - 18 ergs -1 cm - 2 arcsec -2 . Although the instruments used and depths achieved
are heterogeneous, all of the aforementioned observations reached sufficient depth to detect
Lya nebulae with SB Lya > 10 - 17 ergs -1 cm -2 arcsec -2 . Taking all of these observations to¬
gether, we have observed a total of 26 quasars and detected a single giant Lya nebula around
UM287 (41), implying that 1/26 = 0.04 of quasars host giant Lya nebulae. Adopting the
la confidence interval appropriate for the Poisson distribution in the small number regime, we
deduce that the frequency of occurrence of giant nebulae is in the range 3 — 9%.
Finally, an independent narrow band imaging survey of quasars at z ~ 2.7 is being car¬
ried out as part of the Keck Baryonic Structure Survey (KBSS), and the first results on the
distribution of LAEs around quasars were recently published by Trainor et al. (118). Based
on these narrow band imaging observations, (39) published the discovery and additional ob¬
servations of a giant Lya nebula around the quasar HS1549+19 with properties comparable to
that around SDSSJ0841+3921. Trainor et al. quote a typical narrow-band imaging depth of
m NB (3a) 26.7 for point sources, resulting from 5-7hr integrations on Keck LRIS using filters
with FWHM r^j 80A. This corresponds to a (la) SBL ya — 1.5 x 10 18 ergs 1 cm 2 arcsec 2 ,
which although shallower than our own Keck LRIS survey (owing to their wider filters), it
is nevertheless sufficient to detect giant nebulae with S 11 r jya > 10 -17 ergs -1 cm -2 arcsec -2 .
There are 8 objects observed in the Trainor et al. sample, but it is unknown to us whether ad¬
ditional Lya nebulae are present, besides that around HS 1549+19. As such, one can place a
lower limit on the frequency ofgiantLya nebulae to be 1/8 = 0.125, ora la Poisson confidence
interval of 10 — 29%.
To summarize, three independent surveys indicate that the frequency of occurrence of giant
Lya nebulae around quasars is consistent with being ~ 10%, although there are still limitations
due to relatively small samples, heterogeneity of the datasets, and the differences in the search
techniques (narrow-band imaging versus spectroscopy).
8 The Environments of Other Quasars with Giant Lyci Neb¬
ulae and the Impact of Lya Fluorescence
Although quasar clustering indicates that the majority of z ~ 2 quasars reside in moderate
overdensities, we speculate that ~ 10 % trace much more massive structures, and it is this subset
which also exhibits giant Lya nebulae. SDSSJ0841+3921 clearly supports this interpretation,
but there are indeed other examples.
Of course, we also have NB imaging data for UM287, the quasar at z — 2.27 with a gi¬
ant Lya nebula recently discovered by our group (41). However, because UM287 is a hyper-
luminous quasar, we cannot characterize its environment via the number counts of LAEs as we
have for SDSSJ0841+3921. The reason for this is that distribution of LAEs on Mpc scales
from hyper-luminous quasars is likely to be dramatically enhanced by Lya fluorescence (61),
and it is not straightforward to disentangle the fluorescence and clustering enhancements. For
example, consider the hyper-luminous quasar HE0109—3518, the proto-typical case for Lya
fluorescence by dark-galaxies studied by (61). The number of LAEs within R < Ah~ l cMpc
of HE0109—3518 is ~ 8.5 times larger than the cosmic background number expected from the
LAE luminosity function, whereas the quasar-LBG correlation function (a reasonable proxy for
the quasar-LAE clustering) averaged over the volume probed by these narrow-band observa¬
tions would only predict an enhancement of ~ 1.5. One comes to similar conclusions about the
number of LAEs discovered around the hyper-luminous quasars in the KBSS sample of (118).
Progress on understanding the environment of UM287, and its relationship to the giant Lya
nebulae, will require a sample of galaxies selected via broad-band photometry (i.e. LBG selec¬
tion), to avoid contamination from fluorescence. Nevertheless, while fluorescence precludes a
simple interpretation of LAE clustering around UM287, we note that, intriguingly, follow-up
observations also revealed a faint V ~ 22 companion quasar, indicating an overdensity of AGN
similar to SDSSJ0841+3921.
Note however that enhanced LAE clustering around SDSSJ0841+3921 due to Lyo; fluores¬
cence is expected to be negligible, because brightest quasar in the system, namely the f/g quasar
with V = 19.8 is an order of magnitude fainter than the hyper-luminous quasars which exhibit
this fluorescence effect. Radiative transfer simulations of Lya fluorescence around quasars (61)
supports this conclusion. Furthermore, the four AGN in SDSSJ0841+3921 clearly constitute a
dramatic overdensity of AGN, and there is no way for Lya fluorescence to enhance AGN clus¬
tering. Note that the other three AGN do not significantly enhance the ionizing luminosity and
hence cannot boost the fluorescence signal. The fainter Type-1 AGN AGN2 and AGN3 have
V ~ 23 (see Table S2), and thus contribute negligibly to the total ionizing photon production
compared to the f/g quasar. Although we do not detect the continuum from the accretion disk
of the Type-2 AGN AGN1, it is detected in the WISE bands W1 and W2 (see Table S2). At
z ~ 2 these bands probe rest-frame ~ 1 //in. which are likely still dominated by emission from
the reddest part of the accretion disk. The fact that AGN 1 is > 1 magnitude fainter than the f/g
quasar indicates that its UV ionizing continuum (although obscured from our perespective) is
very unlikely to be brighter than that of the f/g quasar.
The KBSS quasar HS1549+19 at z — 2.85, which has been narrow band imaged by Trainor
et al. (118), harbors a giant Lya nebula (39) with SB and properties comparable to that around
SDSSJ0841+3921. Extensive imaging and spectroscopy of HS1549+19 has demonstrated
that it is clearly a dramatic protocluster, with a Mpc-scale overdensity of Lyman break galaxies
(LBGs) a factor of ~ 3.5 times larger than that around typical quasars (10, 38). Although this
quasar is hyper-luminous, Lyct fluorescence is not an issue because the protocluster overdensity
is based on the clustering of broad-band selected LBGs, and not LAEs.
Finally, one may have concerns that the dramatic enhancements of LAEs around HzRGs
and LABs, quantified by our clustering analysis in Supp. 4 and shown in Figs. 2 and S5, could
be partly due to Lya fluorescence. This is a valid concern, because there is compelling evidence
that the HzRGs are obscured quasars (119), and their ionizing photons, although obscured from
our perspective, surely shine in other directions, and are likely to power their giant Lya nebu¬
lae. Similar considerations hold for LABs, although a direct association with obscured AGN has
been more challenging to establish (18,120-123). However, we think that it is unlikely that flu-
orescence plays a significant role in the strong clustering of LAEs around HzRGs and LABs for
two reasons. First, there is no evidence that the unobscured counterparts of HzRGs and LABs
(i.e. the Type-1 quasar shining in other directions) are as bright as the hyper-luminous V < 17
quasars, and hence sufficient to power LAE fluorescence on Mpc scales. Furthermore, even if
the HzRGs and LABs were shining brightly in unobscured directions, the LAE enhancement is
expected to be significantly smaller than for unobscured quasars, because given the cylindrical
volume probed by NB imaging observations, the illuminated volume will be much larger for an
unobscured quasar pointing toward the observer and hence along the major axis of the cylindri¬
cal survey volume, than for an obscured quasar which illuminates a volume randomly oriented
relative to our line-of-sight direction (61).
Thus, although based on only a handful of objects, existing observations support a picture
where ~ 10% of radio-quiet quasars are surrounded by giant Lya nebulae, which acts as a
signpost for the presence of a protocluster. A similar connection between AGN activity, giant
Lya nebulae, and protocluster overdensities has been suggested by past work on HzRGs and
LABs. Our clustering analysis in Supp. 4 (see Figs. 2 and S5) demonstrates that these overden¬
sities are indeed very large and well in excess of the environment of radio-quiet quasars. These
protocluster overdensities are real, and are not the result of Lya fluorescence.
9 Absorption Line Analysis
Given the experimental design of our QPQ survey, we are able to explore characteristics of
the gas surrounding the quadruple AGN in unprecedented detail. Fig. 4 shows strong Hi Lya
absorption coincident with the redshift of the protocluster. We measure a rest-frame equivalent
width W\ = 3.75±0.05 A, offset by a velocity Sv = cAA/A ~ +650 km s” 1 from the f/g quasar
systemic redshift. The absence of very strong damping wings limits the HI column density to
be IVhi < 10 19 7 cm~ 2 . We also present our favored model which has N m = 10 19 2± °' 3 cm -2 ,
suggesting optically thick gas at the Lyman limit. It is possible, in the presence of significant
line-blending, that the Nm value could be much lower, but we consider this scenario very
unlikely, since it would imply a gas metallicity exceeding the solar abundance (see below).
Fig. S7 presents spectra of the b/g quasar centered at the wavelengths of a series of strong,
commonly observed UV resonance-line transitions (including the several transitions shown in
Fig. 4 of the Main text) in the rest-frame of the f/g quasar (based on its near-IR systemic red¬
shift). The gas exhibits especially strong low-ion absorption (Ol 1302, Cll 1334) and moder¬
ately strong intermediate ion absorption (Nil 1083, Sim 1206). The latter absorption occurs
within the Lya forest and could be partially blended with coincident IGM absorption, but we
believe such contamination is modest because the line-profiles track one another as well as the
other low-ion transitions. Surprisingly, there is no statistically significant absorption at the high-
ion transitions of Si +3 or C +3 . Even without performing a detailed photoionization analysis, we
can already conclude that the gas is not highly ionized, which further suggests the gas is opti¬
cally thick to ionizing radiation. One draws a similar conclusion from the relative strengths of
OI 1302 to Sill 1304, where the latter is observed to have larger equivalent width in systems
that are highly ionized (e.g., (124,125)), contrary to our measurements for this system.
For each of the absorption lines, we have estimated rest-frame equivalent widths with sim¬
ple boxcar integration. Ionic column densities were estimated from the metal-line transitions
assuming the linear curve-of-growth approximation and adopting lower limits to positive de¬
tections and upper limits for non-detections (Table S6). Because neutral oxygen 0° undergoes
charge-exchange reactions with atomic hydrogen, the ratio N 0 o /N m is very nearly equal to the
number ratio O/H over a wide range of physical conditions and ionizing radiation fields (25).
Therefore, we estimate a conservative lower limit to the oxygen metallicity [O/H] > —1. To
estimate abundances of the other ions and to constrain physical conditions of the gas, we must
evaluate the ionization state with photoionization models (see the next section).
Our previous studies have shown that cool gas with strong metal-line absorption from lower
ionization states is common in the environment surrounding z ~ 2 quasars (26, 28). In Fig. S8
we show the distribution of equivalent widths for Hi Lya and Cll 1334 from the set of QPQ
pairs with separations less than 200 kpc. The absorption strength for SDSSJ0841+3921 lies
toward the upper end of the distribution but is not extreme.
In Table S6 we include estimates for the abundances of the observed elements assuming
ionization corrections from the photoionization model discussed below and N m = 10 19 2 cm -2 .
The observations require >1/10 solar metallicity and several measurements suggest solar (or
even super-solar) enrichment.
10 Joint Constraints on Protocluster Gas in Absorption and
Our preferred mechanism for powering the giant Lya nebula in SDSSJ0841+3921 is fluo¬
rescent emission from gas in the CGM of the protocluster powered by the radiation from the
f/g quasar, with the gas residing in a population of cool dense clouds which are optically thin
(IVhi < 10 17 ' 5 ) to ionizing radiation. Because our absorption line observations of the b/g quasar
provide constraints on the gas properties in the protocluster, we adopt a novel approach of
modeling the emission and absorption properties of this gas simultaneously, using the Cloudy
photoionization code (126). First, we describe our model for the distribution of the cool clouds,
and introduce simple analytical models for the fluorescent emission to build intuition about the
scaling with gas properties. Then, we consider a sequence of photoionization models of the gas
intercepted by the b/g sightline, which is however undetected in emission. Our absorption line
measurements and photoionization modeling imply that the total cloud column density must be
log 10 i+H = 20.4 ± 0.4. Under the reasonable assumption that the Lyo; emitting gas in the neb¬
ula has the same Am as that detected in absorption, we then show that reproducing the emission
level requires very high gas volume densities tir — 2.0 cm -3 .
10.1 Cool Cloud Model and Emission Estimates
We consider an idealized model whereby the gas in the protocluster halo resides in a popula¬
tion of cool T ~ 10 4 K spherical gas clouds, which have a single uniform hydrogen volume
density ?z H and hydrogen column density iV H . The clouds are uniformly distributed throughout
a spherical halo of radius R, and their aggregate covering factor, defined to be an average over
the sphere, is f c . Four parameters (-R,71 h ,-/Vh, /c) completely describe the distribution of the
gas (see ( 22 ) for further details), for example the total cool gas mass can be written
M c = 3.3 x 10
250 kpc J \ 10 20 - 5 cm -2
— ) m q
0.5 1 J
We now consider fluorescent Lya emission from this population of cool gas clouds powered
by a bright central quasar. There are two regimes where simple expressions can be obtained for
the resulting emission, namely when the gas clouds are either optically thin N m < 10 17 5 cm -2
or optically thick N m > 10 l7 " 5 cm -2 to Lyman continuum photons. For the optically thin case,
the gas clouds are highly ionized, and the average Lya SB from the halo will be
SB Lya = 1.2 x IQ’ 17
1 + z
1.0 cm -3 /\ 10 20 - 5
erg s 1 cm 2 arcsec 2 .
( 12 )
Note that because the gas is already highly-ionized, the SB is independent of the luminosity of
the central ionizing source provided that it is bright enough to maintain the gas as optically thin.
If the clouds are optically thick to ionizing radiation, they will no longer emit Lya propor¬
tional to the volume that they occupy because of self-shielding effects. Instead, a thin highly
ionized skin will develop around each cloud, and nearly all recombinations and resulting Lya
photons will originate in this skin. The cloud will then behave like a ‘mirror’, converting a frac¬
tion 7 / thi ck = 0.66 of the ionizing photons it intercepts into Lyct photons (127). In this regime,
it can be shown (see (22)) that the fluorescent SB, averaged over the halos is
= 8.8 x 10- 17
1 + z
0.5 J V100 kpc
10 3 °-5 ergs -1 Hz
ergs 1 cm 2 arcsec 2 ,
where hv\ jyn is the energy of a Lya photon and <f> (phot s 1
cm 2 )
is the ionizing photon number
To obtain the second equality in eqn (  ), we assume that the quasar spectral energy distribu¬
tion obeys the power-law form L u = L l/LL (i//r' LL ) Q?uv , blueward of i'll and adopt a slope of
ctuv = —1-7 consistent with the measurements of (128). The quasar ionizing luminosity is then
parameterized by L Uhh , the specific luminosity at the Lyman edge. As we describe in detail be¬
low, we estimate that f/g quasar has log 10 = 30.5, which is the value adopted in eqn. ( |T3| ).
In the optically thick regime, we see that the SB now depends on the luminosity of the quasar,
as well as on the covering factor fc, but is independent of the gas properties (i.e. n^andNu).
The smooth morphology of the emission in the Lya nebula implies a covering factor of
fc 2d 0.5. Even if the emitting clouds have angular sizes much smaller than the PSF of our
NB imaging, the morphology of the nebular emission would appear much more clumpy for
lower covering factors. We have explicitly verified this by generating simulated images of
mock Lya emission distributions with the required SB and a range of covering factors fc,
and then convolving them with our NB-imaging PSF (129). In what follows we will always
assume a covering factor of fc = 0.5 for both the gas we detect in emission, as well as that
detected in absorption. This value is well motivated: a) for the emission it is based on the
smooth morphology of the nebula; b) for the absorption it is approximately the covering factor
of optically thick absorbers on R < 200 kpc scales in the quasar CGM (27).
The average SB of the giant Lya nebula at R ~ 100 kpc is SB LyQ! = 1.3 x 10 -17
ergs -1 cm -2 arcsec -2 , which we compute via the emission from the area within our 2 a SB
contour which intersects the the annulus R = [75,125] kpc. This is several times lower than the
expectation from eqn. (  ) for the optically thick regime. Thus our simple analytical scaling
relations already appear to favor a situation where the clouds are optically thin. We return to
this comparison below using photoionization models.
10.2 Details of Cloudy Modeling
The expressions in eqns. < [I2| ) and ( fl3| ) provide reasonable estimates for the fluorescent SB, but
they neglect other sources of Lya emission such as collisional excitation (i.e. cooling radiation),
and pumped or ‘scattered’ Lya- photons arising from the diffuse continuum produced by the gas
itsel. Note that we do not attempt to model scattering of Lya photons produced by central f/g
quasar itself. Radiative transfer simulations of radiation from a hyper-luminous quasar through
a simulated gas distribution (41), have shown that the contribution of scattered Lya line photons
from the quasar do not contribute significantly to the Lya surface brightness of the nebula.
Given that this calculation was for a hyper-luminous quasar a factor of ~ 10 times brighter than
the f/g quasar, scattered line photons from the quasar are expected to contribute negligibly to
the nebular emission., and they also do not treat the temperature dependence of the emission.
To fully account for these effects, we construct photoionization models for a slab of gas being
illuminated by a quasar using the Cloudy software package (126). Cloudy predicts the emergent
spectrum from the gas, allowing us to compute its SB, which we then dilute by the covering
factor fc = 0.5.
We model the quasar spectral-energy distribution (SED) using a composite quasar spectrum
which has been corrected for IGM absorption (128). This IGM corrected composite is important
because it allows us to relate the g-band magnitude of the f/g quasar to the specific luminosity
at the Lyman limit L IJ] ] . For energies greater than one Rydberg, we assume a power law form
L v = L UhL (u / u LL ) avv and adopt a slope of «uv = —1-7, consistent with the measurements
of (128). We determine the normalization L Jy| l by integrating the Lusso et al. (128) composite
spectrum against the SDSS filter curve, and choosing the amplitude to give the correct g-band
magnitude of the f/g quasar (see Table S2), which gives a value of log 10 L VLL = 30.5. We
extend this UV power law to an energy of 30 Rydberg, at which point a slightly different power
law is chosen a = —1.65, such that we obtain the correct value for the specific luminosity at
2 keV L u ( 2keV) implied by measurements of a ox , defined to be L v (2 keV)/Lj,(2500 A) =
(^ 2 kev/^ o) Q ° x . We adopt the value a ox = —1-5 measured by (130) for SDSS quasars
of comparable luminosities. An X-ray slope of a x = — 1, which is flat in vf v is adopted
in the interval of 2-100kev, and above 100 keV, we adopt a hard X-ray slope of o H x = —2.
For the rest-frame optical to mid-IR part of the SED, we splice together the composite spectra
of (128), (84), and (131). These assumptions about the SED are essentially the standard ones
used in photoionization modeling of AGN (e.g. (132)). We also include the contribution from
the radiation field of the ambient UV background radiation field using the spectrum of (133),
but it is completely sub-dominant relative to the f/g quasar at the distances that we consider.
These Cloudy models require as an input the ionizing photon number flux <F(r) (see
eqn. (fl~4|)) at the surface of the cloud. Because this varies with radius as r” 2 , in principle mod¬
eling the halo emission would require integrating over a distinct photoionization model at each
radius. For simplicity, we restrict attention to just a few radii. For our model of the optically
thick absorber, we set r = 250 kpc, which we believe is the likely location of the absorbing
gas for several reasons. First, the nebula itself has an end-to-end size of ~ 500 kpc roughly
centered on the f/g quasar, such that r = 250 kpc is a reasonable radius for the absorber given
the measured impact parameter R± = 176 kpc of the b/g quasar sightline. Second, given this
impact parameter, the characteristic velocity separation of the absorbing gas from systematic
Umax = 600 kms -1 implied by the metal-line absorption (see Fig. S7), we can compute the con¬
ditional probability distribution P(< r\R±, u max ) that the distance between the quasar and the
absorber is less than r, which is determined by the quasar-LLS correlation function which we
have measured (24, 27). For the problem at hand, we find that P(< r\R±,v ma , x ) = 0.43 for
r = 250 kpc, again indicating that this is a reasonable distance (see a similar discussion in § 4.1
of (25)). According to eqn. 0. the SB arising from such optically thick gas scales with <f>,
and hence will decrease as r~ 2 . For the gas producing the Lya emission, we model three radii
r = [50,100, 250] kpc. We will argue that the emitting gas in the nebula is optically thin, and in
this regime the SB is independent of $ (see eqn. (121). Thus considering a few radii to describe
the emission is a reasonable approximation provided the gas remains optically thin at the radii
The Cloudy model of the cool cloud is completely specified by the value of <f>, the SED
shape, the metallicity of the gas Z, the hydrogen volume density n H , and a ‘stopping criterion’,
which sets the total column density of the cloud. Because we assume the f/g quasar and UVB
are the only sources of ionizing photons, <I> and the SED shapes are fixed. For the metallicity,
we adopt a value of log (Z/Z Q ) = —0.5 consistent with our lower-limit (log (Z/Z Q ) > —1) on
the metallicity of the optically thick absorber detected in the b/g quasar sightline (see Supp. 9).
We have confirmed that the results discussed below are insensitive to the assumed gas metal¬
licity. We have also verified that our results are insensitive to the exact value of the UV slope
a = —1.7, within the 1 a error measured by Lusso et al. a = —1.7 ± 0.6. Our models thus
constitute a one-dimensional sequence parameterized by the volume density of the clouds nn.
Photoionization models are known to be self-similar in the ionization parameter U = T /riuc,
which is the ratio of the number density of ionizing photons to hydrogen atoms. Since $ has
been fixed, our sequence of models is thus also a sequence in U, spanning the relevant range of
ionization conditions. Our stopping criterion depends on what is known about the gas that we
are attempting to model, which we elaborate on next.
10.3 Photoionization Model of the Absorbing Gas
Here we compare photoionization models to the properties of the optically thick absorber in
the spectrum of the b/g quasar, lying at at an impact parameter of R± = 176 kpc from the f/g
quasar. We follow the standard approach (25) and focus our analysis on multiple ionization
states of individual elements to avoid uncertainties related to intrinsic abundances. Specifically,
our column density measurements in Table S 6 yield the following limits on ionic ratios:
C+/c 3 + > 0.8, Si + / Si 3+ > 1.3, and Si + /Si 2+ < 0.7. In addition, we measured the neutral
column density to be Am = 10 19 ' 2 cm -2 , which we use as the stopping criterion of the Cloudy
model. In other words, Cloudy will continue to add slabs of material until this neutral column
is reached, guaranteeing that we satisfy this constraint. Finally, we do not detect Lyo emission
from the vicinity of the background sightline, allowing us to set an upper limit of SB L y Q < 5.1 x
10 -18 ergs -1 cm -2 arcsec -2 from the cool gas. This SB limit was computed via simulations
whereby a fake circular source with radius of 3" was placed at the location of the b/g quasar.
This exercise showed that only sources with SB LyQ > 3 x SB lcr would be clearly detectable,
where SB lcr = 1.7x 10 -18 ergs -1 cm -2 arcsec -2 is our 1 a SB limit for 1.0 sq. arcsec apertures.
Note that the individual clouds intercepted by our sightline may have extremely small sizes
~ 10 pc making it impossible to angularly resolve emission from a single cloud. However, we
have argued that these clouds have a large covering factor of / c = 0.5 such that we would be
able to detect their aggregate emission near the location of the b/g sightline (see Fig. 1 in (22)
for more details).
Our sequence of photoionization models is compared to the observational constraints in
Fig. S9. The three different ionic ratios considered paint a consistent picture for the ioniza¬
tion state, and we constrain the ionization parameter to be log 10 U = —2.7 ± 0.3, where we
conservatively include an error of 0.3 dex to account for both statistical and systematic uncer¬
tainties in the photoionization modeling. Accordingly, the neutral fraction is constrained to be
log 10 xhi = —1-2 ± 0.3, and the total hydrogen column log 10 A H = 20.4 ± 0.4, where the er¬
ror on Ah includes the quadrature sum of modeling uncertainty in the neutral fraction, and the
uncertainty in our measurement of Ahi (see Supp. 9). Plugging this determination of the total
hydrogen column density into eqn. ([IT]) we obtain a total cool gas mass of M c = 2.4 x 10 11 M 0
within a spherical halo R = 250 kpc, or the range 1.0 x 10 11 M 0 < M c < 6.5 x 10 11 M 0 im¬
plied by the la measurement and modeling error on A r tI . Note that given the self-similarity of
the photoionization models with ionization parameter U, these results are essentially indepen¬
dent of the value of $ that we have assumed in our models, and hence the distance of the gas
from the quasar. Lastly, we have calculated the ionization corrections to the ionic column den¬
sities measured for the absorbing gas and provide limits on the elemental abundances implied
by each metal absorption line measurement in Table S 6 .
According to the photoionization models, if indeed the absorbing gas lies at a distance
r = 250 kpc and is illuminated by the f/g quasar, the Lya emission would be SB LyQ ~
3 — 4 x 10 -18 erg s -1 cm -2 arcsec -2 , which lies just below our detection limit SBL yQ! = 5.1 x
10 -18 ergs -1 cm -2 arcsec -2 , and is hence consistent with our non-detection of extended emis¬
sion near the b/g quasar sightline. Indeed, because the gas is optically thick, and essentially
absorbs all ionizing photons that impinge on it, the SB LyQ! is nearly independent of the volume
density, and as expected from eqn. (  ), depends primarily on $ and our assumed covering fac¬
tor fc = 0.5. The weak dependence on n H arises because, at the lowest ionization parameters
and hence highest neutral fractions, as much as ~ 20% of the Lyo is produced by collisional
excitation (green curve labeled collisions in Fig. S9). This is because collisional excitation of
neutral hydrogen requires a large neutral fraction, and it is also exponentially sensitive to tem¬
perature. The low ri]j < 10 - 2 cm 3 and high ionization parameter log 10 U > — 1 models are
highly ionized xhi < 10 -3 , and have higher temperatures T > 17, 000 K, both of which tend
to suppress collisional excitation. As a result of the very weak dependence of SB LyQ on ri\\ in
the optically thick regime, our non-detection of emission does not allow us to put interesting
constraints on the volume density n H for the absorbing gas.
Note that the fact that our photoionization model predicts Lyo- emission from the absorbing
gas just below our detection limit results from our assumed distance of r = 250 kpc, since
for optically thick gas SB LyQ oc </> cx r -2 (see eqn. 13). Had we assumed a distance equal
to the impact parameter R± = 176 kpc the predicted emission would be a factor of two higher
SB LyQ ~ 10 -17 ergs -1 cm -2 arcsec -2 , and hence detectable. As such it is informative to briefly
discuss all possible scenarios that would result in a non-detection of Lyct from the absorbing
gas. These are: 1) the gas is illuminated by the f/g quasar, but it lies at a distance r > 250 kpc
such that its fluorescent Lya emisssion is below our detection threshold; 2) the gas is illuminated
by the quasar and is at a distance comparable to the impact parameter, but the covering factor of
optically thick gas is much lower than we have assumed fc = 0.5; 3) the absorbing gas is not
illuminated by the f/g quasar but is instead shadowed, perhaps because from its vantage point the
quasar is extincted by the same obscuring medium invoked in AGN unification models (67,68).
Unfortuately, our current observations cannot distinguish between these three scenarios,
however we favor scenario 1), that the gas lies at a distance ~ 250 kpc from the quasar, as
shown in Fig. S9. While it is possible that the covering factor of optically thick gas is indeed
lower than the f c = 0.5 that we have assumed, we consider this unlikely for two reasons.
First, this would be in conflict with the statistical properties of optically thick absorption for the
aggregate population of quasars (27). Although, SDSSJ0841+3921 is atypical given that it is
a protocluster and hosts a giant nebula, we already argued that the smooth morphology of the
giant nebula also argues for a high covering factor f c > 0.5 of the emitting gas, to prevent the
emission from appearing too clumpy. Under the reasonable assumption that the absorbing and
emitting gas have the same origin, one expects their covering factors to be comparable.
At first glance, the possibility that the absorbing gas is shadowed from the quasar radiation
seems quite plausible. Indeed, we have previously argued in the QPQ series that the anisotropic
clustering of LLSs around quasars implies that the optically thick gas in the quasar CGM is
often shadowed (22-24). However, we robustly determine the ionization parameter of the ab¬
sorbing gas to be log 10 U = —2.7 ± 0.3. If obscuration shadows the gas from the quasar, then
the radiation field would be expected to be the UV background, which given our ionization
parameter, would then imply a gas density riu ~ 10 - 3 cm~ 3 . In the next section, we will see
that the detection of bright Lya emission implies that the density of the emitting gas is three
orders of magnitude higher n H ~ 1 cm -3 . As indicated in the middle panel of Fig. S9, if the
gas is illuminated by the quasar at a distance r = 250 kpc, then gas density implied by our ion¬
ization parameter is then % — 0.7 cm -3 , which is then consistent with the value deduced from
modeling the emission. Thus again invoking the very reasonable assumption that the absorbing
and emitting gas have the same origin, the fact that we infer a comparable riu for the absorb¬
ing gas as required for the emitting gas, seems to favor a scenario where the absorbing gas is
illuminated by the f/g quasar, but lies at a distance r > 250 kpc, such that the fluorescent Lya
is below our threshold. Future deeper Lya observations could test this hypothesis if they can
reach SE>Ly a — 10 - 18 ergs _ 1 cm _ 2 arcsec -2 levels, such that one could then detect emission
from the absorbing gas out to a distance as far as r ~ 500 kpc.
10.4 Photoionization Model of the Emitting Gas
Here we consider photoionization models for the emitting gas under the assumption that it has
a comparable column density jV H to the gas we detect in absorption. In modeling the emission
we consider gas at r = 100 kpc which is the characteristic scale of the emission nebula. We
assume that the total hydrogen column density does not vary strongly between these distances.
This hypothesis is supported by a photoionization modeling study of a large sample of quasar
CGM absorbers, which indicates that on average the total column density varies weakly (if
at all) with impact parameter on such scales. We run the same Cloudy photoionization models
as described previously, but with the modification that the stopping criterion is now chosen to
be that the total hydrogen column A'n = 10 20 4 cm -2 , deduced in the previous section.
In Fig. S10 we show the results of our sequence of photoionization models for the Lya
emitting gas. The black curve illustrates the dependence of SB LyQ , on volume density n H for
our fiducial choice of the radiation field intensity, where the gas lies at a distance of r = 100 kpc.
The upper x-axis indicates the corresponding ionization parameter U at this distance. The red
and blue curves are for stronger and weaker radiation fields, corresponding to r = 50 kpc
and r = 250 kpc, respectively. For low values of «h the gas remains optically thin A r ni <
10 17 ' 5 cm -2 , and the SB Lya exhibits a nearly linear dependence on the volume density n H and
is nearly independent of the ionizing radiation intensity, as predicted by eqn. 0- The gas
becomes self-shielding and optically thick once N m > 10 1 ' 5 cm -2 , at which point the SB LyQ
plateaus at a single value depending only on <1> (since we have fixed the covering factor fc =
0.5) and hence on distance r, as predicted by eqn. ( fT3| ).
These models clearly indicate that, given the column density of cool material Nu = 10 20 4
cm -2 , the volume density must be % — 2.0 cm” 3 , in order to reproduce the observed average
emission level at 100 kpc of SB Lya = 1.3 x 10” 17 ergs -1 cm” 2 arcsec” 2 (indicated by the
horizontal dashed line in Fig. S10). Provided that the emitting gas is optically thin N m < 10 17 5 ,
this conclusion is essentially indepdenent of the radiation field and hence the assumed distance
of the gas from the f/g quasar. In this regard, the relatively flat radial SBi, ya profile observed for
the giant nebula is readily explained. It naturally arises because order of magnitude variations
in the radiation field, corresponding to gas at distances ranging from r = 50 — 200 kpc, produce
nearly the same level of Lyo emission, provided that IVh does not vary significantly across the
nebula, as we have assumed.
A corollary of our absorption-line constraint on Nn — 10 2 ° 4 cm” 2 and our emission con¬
straint on tt-h — 2.0 cm” 3 is that the implied absorption path length through the cool clouds is
^cloud = N K /n K rs-/ 40 pc, implying that the emitting gas is distributed in extremely compact
A v (km/s)
-4000 -2000 0 2000 4000
-4000 -2000 0 2000 4000
A v (km/s)
A v (km/s)
-4000 -2000 0 2000 4000
-4000 -2000 0 2000 4000
A v (km/s)
Fig. SI: Two dimensional spectrum plotted as a x-map for the three slits in Fig. 3. The lower and
middle rows of each image show Xsky (sky-subtracted only) and Xsky+PSF (sky and PSF subtracted),
respectively as defined in § |T] (see ( 22 ) for additional details). The upper row shows the smoothed map
Xsmth? (sky and PSF subtracted, then smoothed), which is helpful for identifying extended emission. A
symmetric Gaussian kernel (same spatial and spectral widths in pixels) was adopted, with dispersion
<7 sm th = 100 km s ' 1 (FWF1M = 235 km s -1 ), and ler spatial smoothing of 0.65" (FWFLM^l.5")- In the
absence of extended emission, the distribution of pixel values in the the sky and PSF subtracted images
should be Gaussian with unit variance. The lower and middle rows are displayed with a linear stretch
ranging from — 6a to 6 cr indicated by the color-bar at middle right. The upper row has a stretch from
—6a to 16<7 is adopted, and a 7 value chosen to decrease the range of color at the highest significance, as
indicated by the color-bar at upper right. The horizontal dashed lines indicates the spectral traces for the
f/g and b/g quasar, and the vertical dashed line indicates An = 0, corresponding to the systemic redshift
of the f/g quasar. Note that the b/g quasar spectrum in Slit3 (right) has much lower S/N ratio than in
Slit2 (center) because of much poorer seeing.
Fig. S2: Optical spectra of the four AGN at the same redshift. From top to bottom: f/g quasar,
AGN1, AGN2, and AGN3. The gaps in the coverage of the spectra of AGN1 and AGN2 are due to the
LRIS configuration chosen to cover both Lyo, C IV, and He II. The spectrum of AGN3 does not cover
the Lya: line due to the limited blue sensitivity of DEIMOS. Note the different scale on the y-axis to
better show the different emission lines, which are marked with vertical dotted lines and labeled.
Fig. S3: Comparison of the Line-ratios of AGN1 to other Type-2 AGN. The line-ratios of AGN1
(red square) are compared to other Type-2 AGN compiled from the literature in the C IV /He II vs C I v/
C III] plane (left) as well as the C III] /He II vs C III] /C II] plane (right). The circles are individual mea¬
surements of HzRGs (black) from the compilation of (80), and narrow-line X-ray sources (cyan) and
Seyfert-2s (blue) from the compilation of (81). Triangles indicate measurements from the composite
spectra of HzRGs (orange) from (70), the composite Type-2 AGN spectra (purple) from (82), who split
their population into two samples above and below Lya EW of 63A, and a composite spectrum of mid-IR
selected Type-2 AGN (green) from (83). The stars represent measurements from composite spectra of
Type-1 quasars, based on the analysis of (magenta) (84), (red-magenta) (85), and (blue-magenta) (86).
The line-ratios for AGN1 lie in the locus of points spanned by other Type-2 objects in the literature.
Fig. S4: Confidence regions on correlation function parameters from maximum likeklihood clus¬
tering analysis. The 1 ,2, and 3a confidence regions are illustrated by the black, blue, and red contours,
respectively. The maximum likelihood value is ro = 29.3 /i -1 Mpc and 7 = 1.5, which is indicated
by the white cross. Due to the low signal-to-noise ratio of the clustering measurement, and the form
adopted for the correlation function £(r) = (r/ro ) -7 (such that amplitude and slope are not indepen¬
dent parameters), a strong degeneracy exists between ro and 7 , such that neither is individually well
Fig. S5: Binned correlation function of LAEs around HzRGs and LABs. The dimensionless corre¬
lation function 7 defined by eqn. (|9|) is indicated by the data points, with la error bars, determined by
taking the 16th and 84th percentiles of the distribution of x resulting from our bootstrap procedure. The
solid red line indicates the best-fit protocluster-LAE cross correlation function (ro = 29.3 h~ 1 Mpc and
7 = 1.5) determined by our unbinned maximum likelihood procedure. The gray shaded region indicates
the la error on our cross-correlation function fit, implied by the erorrs on the parameters ro and 7 . This
region is determined by evaluating 7 for each sample of the boostrap error distributions of ro and 7 , and
taking the 16th to 84th percentile confidence region of 7 . Our best fit correlation function (red line), and
its quoted errors (gray shaded) clearly provides an acceptable representation of the binned data points
and their respective errors.
Fig. S6: Near-IR spectrum of the f/g quasar SDSSJ084158.47+392121.0 (top) and of the b/g quasar
BOSSJ084159.26+392140.0 (bottom) obtained with the NIRI spectrometer on the Gemini-North
telescope. Each quasar shows the detection of strong [O III] emission and modest H/3 emission. The
[O iii]A5008A doublet lines were centroided, giving observed wavelength centroids of A 0 bs = 1-5231 ±
0.0001/mi and 1.60954=0.00006/rm for the f/g and b/g quasars respectively. This gives systemic redshifts
Zf g = 2.0418 for the f/g quasar and 2.21370 for the b/g quasar. The estimated uncertainty on this
redshift measurements is approximately 50 km s , taking into account both the error in our wavelength
calibration and the estimated uncertainty when one uses [O III] to assess the systemic redshift of a quasar.
This uncertainty on the measurement of the redshift is far lower than the width of our custom narrow-
band filter (~ 2600km s -1 ).
Fig. S7: One-dimensional spectra of the b/g quasar centered at the wavelengths of a series of
commonly observed UV-line transitions in the rest frame of the f/g quasar. The gas exhibits espe¬
cially strong low-ion absorption (OlA1302, SillA1304, CllA1334) and moderately strong intermediate
ion absorption (NIIA1083, Si IIIA1206). The latter absorption occurs within the Lya forest and could be
partially blended with coincident IGM absorption, but such contamination should be modest given that
the line-profiles track one another and also the low-ion transitions. In contrast, there is no statistically
significant absorption corresponding to the SilvA1393 and C IVA1548 transitions indicating the gas is
not highly ionized. See Table S6 for a summary of the properties of these metal absorption lines.
£ °- 10
W Lya (Ang)
■ ■ !
Fig. S8: Equivalent widths of cool gas absorption in the QPQ dataset. Measurements of the equiv¬
alent widths of the Cll 1334 transition against the measurement for Hi Lyct. The measurements are
drawn from our QPQ survey (27, 28 ) and the pair that led to the serendipitous discovery of the quadruple
AGN system is marked separately. The gas probed by the background quasar (at an impact parameter
of 176 kpc) shows strong low-ion absorption characterstic of the full population. This indicates a highly
enriched and optically thick gas.
U (ion param)
IO " 3
IO" 2 3
io 1 a
IO " 3 IO " 2 io 1 io° io 1 io 2
n H (cm" 3 )
Fig. S9: Photoionization model of the optically thick absorber in b/g quasar spectrum. The gas
is assumed to be at r = 250 kpc and exposed to the radiation field of the f/g quasar. The models are
a one-dimensional sequence parameterized by the hydrogen volume density ng on the x-axis, which
given that the radiation field is fixed, is also a proxy for the ionization parameter U shown on the upper
x-axis. Upper panel: Variation of the total hydrogen column Am with nn/U, given that we require the
model to reproduce our measured neutral column density Ami = 10 19 ' 2 . The y-axis on the right shows
the corresponding neutral fraction. Our measured log 10 U (middle panel) in turn implies log 10 Am =
20.4 ± 0.4, indicated by the vertical dotted line. Middle Panel: Variation of three different ionic ratios
(Si + /Si 2+ in red; Si + /Si 3+ in green; C + /C 3+ in blue) with nn/U. Our measured lower/upper limits on
these ratios from Table S6 are illustrated by the colored arrows. We obtain a consistent photoioinization
solution implying log 10 U = — 2.7T0.3, indicated by the vertical dotted line. Lower Panel: Predicted SB
of extended Lya emission from the absorbing gas as a function of nn/U. The red and green curves show
the respective contributions from recombinations and collisional excitation, and the black curve shows
the total Lya SB. The horizontal dashed line with arrows indicates our detection limit for an extended
source near the b/g quasar sightline SBL yQ < 5.1 x 10 -18 ergs -1 cm -2 arcsec -2 . If the absorbing gas
is at r = 250 kpc as we have assumed, our models are marginally consistent with a non-detection of
extended Lya emission at that location. If the gas lies at larger distances r > 250 kpc, the predicted SB
will decreaes as r
Fig. S10: Photoionization model of the Ly a emitting nebular gas. The models assume a total hy-
, deduced from our photoionization modeling of the optically
drogen column density N-r = 10
thick absorber in the b/g sightline (see Fig. S9and § 10.3). The gas is assumed to be exposed to the
radiation field of the f/g quasar at distances of r = 50, 100, and 250 kpc, indicated by the red, black,
and blue curves, respectively. The models arc parameterized by the hydrogen volume density n h on
the x-axis. The ionization parameter U corresponding to the gas at r = 100 kpc (black curves) is de¬
noted on the upper x-axis. Upper panel: Variation of the neutral hydrogen column IVhi with «h. The
y-axis on the right shows the corresponding neutral fraction. Lower panel: Variation of the Lya SB
with «h- For low ng the gas remains optically thin JVhi < 10 175 cm -2 , and the SB SBL yQ varies
linearly with «h and is nearly independent of the ionizing radiation intensity. At higher nn the gas
begins self-shielding A r ni > io 17 - 5 cm 2 , and the SB plateaus at a single value depending only on <f>
and hence on distance r. The models indicate that a value of n h — 2.0 cm -3 , indicated by the ver¬
tical dotted line, is required in order to reproduce the observed average emission level at 100 kpc of
SBL yQ , = 1.3 x 10 -17 ergs -1 cm -2 arcsec -2 , indicated by the horizontal dashed line in Fig. S10. This
conclusion is essentially independent of the radiation field and hence the assumed distance of the gas
from the f/g quasar.
(kms -1 )
08 41 58.47
+39 21 21.0
08 41 58.24
+39 21 29.1
08 41 58.66
+39 21 14.7
08 41 59.10
+39 21 04.6
08 41 59.25
+39 21 40.0
Table SI: AGN properties. The coordinates, redshift z, redshift error a z , angular 9 and transverse separa¬
tion Rj_ from the f/g quasar arc listed for the four AGN in the protocluster, as well as the b/g quasar.
J 'C 00 O h ^
j<n d h ^
> —< — oo o
; © — o< r^
: — cn
m i"0 \o m
s, o o o o o
k o o o o ©
-H -H -H -H -H
11 « ^ M ON «
On -< o os
D» on cn cn cn on
^ ^ « in no no
r o on cn cn cn
-|-| 00 06 °o 06
J 8 A A A A
O 10 in «n O
-H <N <N ri -H
_u os — — — o
3 ^ A A A ®
S 2 2
!> o o dd — o
b ~H -H r-C r~C -H
-U 1^ NO - - NO
n 'P 00 A A
^ 2 12 -
H m 00 no >n
•>, O O CN ©
> © d © in o
b -H -H -H 06 -H
-H in on on - o
JJ GO 00 — A
^ n no 00 no
N • NO NO NO _•
* -Hcd r4-H-H
N n 00 on
A A " ^
CD o <D
cu o .3
I & |
5 • 1-1
*3 m co
c3 O' s
k o ^ ® © ©
-fl on on “H “H
-H 00 (n m m -
« ^ A A ^ ^
o,P cn o q
k O O rn O O
_A -H -H ^ -H -H
s p __ _ q
K O _H — — O
11 ~H cd cd on ~H
3 ^ A A A9
iu cn — cn on co
c O' z Z Z O'
a goO O O so
Z 2a < < < 3
O 05 ,J*j
2 H q
s- ifl M
JT es 7.
Line FWHM Flux W A
(kms -1 ) (10~ 17 ergs -1 cm -2 ) (A)
Lya. 1575 ± 36 76.5 ± 1.5 > 118.2
NV.. — —0.4 ± 1.4 —
CIV.. — 1.3 ±0.8 —
Hell. 1154 ±276 2.4 ± 0.4 >15.8
CIII]. 1434 ±171 2.6 ± 0.3 22.6 ±2.2
CII]. — 1.0 ±0.2 11.2 ±2.6
Table S3: Emision line measurements for the Type-2 quasar AGN1. Line fluxes were measured from the
one-dimensional spectrum of the AGN, and errors determined from the la noise vector. For lines de¬
tected at S/N > 3, rest-frame equivalent widths W\ are computed. If the continuum lies above the la
noise vector these are listed as values, otherwise we quote lower limits. The FWHM for emission lines
for which our discovery spectra had sufficient S/N for a measurement (Lya, Hell, and Cm]) are also
(10 -17 ergs -1 cm -2 )
(10 -19 ergs — 1 cm — 2 A — 1 )
Properties of the candidate LAEs around the f/g quasar. Here / Lya is
the Lya line flux in
erg s 1 cm 2 , f\ is the flux density in
the continuum in
units of 10 19 erg s 1 cm
- 2 A" 1
estimated from the
V-bantl, Wi, ya is the rest-frame Lya equivalent width, and R± the proper projected distance from the f/g quasar in
kpc. LAEs are selected using the procedure described in § 4.1, namely they are required to have Hj, yo > 20A and
fhya > 5.4 x 10 -18 erg s -1 . Objects which were used to measure the overdensity profile around the f/g quasar in
Fig. 2 are denoted by a * in the column Tag, where the name of the source is also indicated. These sources had to
satisfy a brighter flux limit, namely fhya > 4.0 x 10“ 11 erg s -1 cm -2 , and lie within 2.08' of the f/g quasar.
Filter FWHM (A)
EW limit (A)
20 51 03.5
-27 03 04.1
11 40 48.2
-26 29 09.5
14 34 11.0
+33 17 32.6
00 54 29.8
-23 51 31.1
09 45 32.7
-24 28 49.7
22 17 26.0
+00 12 37.7
03 18 12.0
-25 35 10.8
20 09 48.1
-30 40 07.4
Table S5: Properties of objects analyzed in our protocluster-LAE correlation analysis. The nature of the
source (HzRG or LAB), coordinates, redshift z, the FWHM of the NB filter used (A), the limiting EW
of the LAE search (Wij m i t ), and the reference for the observations arc indicated. References used arc:
(1) Venemans et al. 2007 (73); (2) Prescott et al. 2008 (79); (3) Nestor et al. 2011 (106).
(km s -1 )
3750.0 ± 50.0
672.9 ± 14.0
335.0 ± 110.0
Table S6: Absorption line measurements. The columns indicate the Ion, rest wavelength A, oscillator
strength log /, velocity interval n; nt used for the calculation of the rest equivalent widths, rest-equivalent
width W\, ionic column density TV, and elemental abundances [X/H] implied by each metal absorption
line measurement. The column densities TV reported assume the linear curve-of-growth approximation.
The metal abundances [X/H] arc computed using ionization corrections from our favored Cloudy model
presented in Section S10, and arc reported according to the standard convention of logarithmic relative
to Solar (i.e. 0 = solar abundance).
References and Notes
1. J. Magorrian, et al, AJ 115, 2285 (1998).
2. L. Ferrarese, ApJ 578, 90 (2002).
3. J. N. Bahcall, S. Kirhakos, D. H. Saxe, D. P. Schneider, ApJ 479, 642 (1997).
4. M. G. Haehnelt, M. J. Rees, MNRAS 263, 168 (1993).
5. S. Djorgovski, The Space Distribution of Quasars, D. Crampton, ed. (1991), vol. 21 of Astronomical Society
of the Pacific Conference Series, pp. 349-353.
6. J. F. Hennawi, et al.,AJ 131, 1 (2006).
7. S. G. Djorgovski, et al., ApJL 662, LI (2007).
8. E. P. Farina, C. Montuori, R. Decarli, M. Fumagalli, MNRAS 431, 1019 (2013).
9. C. Lacey, S. Cole, MNRAS 262, 627 (1993).
10. R. F. Trainor, C. C. Steidel, ApJ 752, 39 (2012).
11. M. White, et al., MNRAS 424, 933 (2012).
12. N. Fanidakis, A. V. Maccio, C. M. Baugh, C. G. Lacey, C. S. Frenk, MNRAS 436, 315 (2013).
13. B. P. Venemans, et al„ A & A 461, 823 (2007).
14. P. J. Francis, et al„ ApJ 457, 490 (1996).
15. C. C. Steidel, et al, ApJ 532, 170 (2000).
16. A. Dey, et al, ApJ 629, 654 (2005).
17. R. A. Overzier, et al.,ApJ 771, 89 (2013).
18. M. K. M. Prescott, I. Momcheva, G. B. Brammer, J. P. U. Fynbo, P. Mpller, ArXiv e-prints (2015).
19. M. K. M. Prescott, N. Kashikawa, A. Dey, Y. Matsuda, ApJL 678, L77 (2008).
20. Y. Yang, A. Zabludoff, C. Tremonti, D. Eisenstein, R. Dave, ApJ 693, 1579 (2009).
21. Y.-K. Chiang, R. Overzier, K. Gebhardt, ApJ 779, 127 (2013).
22. J. F. Hennawi, J. X. Prochaska, ApJ 766, 58 (2013).
23. J. F. Hennawi, et al„ApJ 651, 61 (2006).
24. J. F. Hennawi, J. X. Prochaska, ApJ 655, 735 (2007).
25. J. X. Prochaska, J. F. Hennawi, ApJ 690, 1558 (2009).
26. J. X. Prochaska, J. F. Hennawi, R. A. Simcoe, ApJL 762, L19 (2013).
27. J. X. Prochaska, et al, ApJ 776, 136 (2013).
28. J. X. Prochaska, M. W. Lau, J. F. Hennawi, ApJ 796, 140 (2014).
29. See further information in supplementary online material § S2 on Science Online.
30. I. Kayo, M. Oguri, MNRAS 424, 1363 (2012).
31. See further information in supplementary online material § S3 on Science Online.
32. R. Ciardullo, et al, ApJ 744, 110 (2012).
33. See further information in supplementary online material § S4 on Science Online.
34. See further information in supplementary online material § S6 on Science Online.
35. See further information in supplementary online material § S5 on Science Online.
36. B. D. Lehmer, et al., ApJ 691, 687 (2009).
37. See further information in supplementary online material § S7 on Science Online.
38. R. E. Mostardi, et al., ApJ 779, 65 (2013).
39. D. C. Martin, et al., ApJ 786, 106 (2014).
40. See further information in supplementary online material § S8 on Science Online.
41. S. Cantalupo, F. Arrigoni-Battaia, J. X. Prochaska, J. F. Hennawi, P. Madau, Nature 506, 63 (2014).
42. M. Fumagalli, et al., ApJ 780, 74 (2014).
43. C.-A. Faucher-Giguere, et al., ArXiv e-prints (2014).
44. S. Lopez, et al., ApJ 679, 1144 (2008).
45. See further information in supplementary online material § S9 on Science Online.
46. See further information in supplementary online material § S10 on Science Online.
47. N. H. M. Crighton, et al., ArXiv e-prints (2014).
48. A. H. Mailer, J. S. Bullock, MNRAS 355, 694 (2004).
49. H. J. Mo, J. Miralda-Escude, ApJ 469, 589 (1996).
50. T. M. Heckman, S. A. Baum, W. J. M. van Breugel, P. McCarthy, ApJ 338, 48 (1989).
51. M. McDonald, S. Veilleux, D. S. N. Rupke, R. Mushotzky, ApJ 721, 1262 (2010).
52. K. N. Abazajian, et al.ApJS 182, 543 (2009).
53. J. F. Hennawi, et al.ApJ 719, 1672 (2010).
54. C. P. Ahn, et al., ApJS 203, 21 (2012).
55. J. B. Oke, et al., PASP 107, 375 (1995).
56. E. Bertin, S. Arnouts, A&AS 117, 393 (1996).
57. E. Bertin, Astronomical Data Analysis Software and Systems XV, C. Gabriel, C. Arviset, D. Ponz, S. Enrique,
eds. (2006), vol. 351 of Astronomical Society of the Pacific Conference Series, p. 112.
58. E. Bertin, et al.. Astronomical Data Analysis Software and Systems XI, D. A. Bohlender, D. Durand, T. H.
Handley, eds. (2002), vol. 281 of Astronomical Society of the Pacific Conference Series, p. 228.
59. J. B. Oke, AJ 99, 1621 (1990).
60. A. U. Landolt, AJ 104, 372 (1992).
61. S. Cantalupo, S. J. Lilly, M. G. Haehnelt, MNRAS 425, 1992 (2012).
62. S. M. Faber, et al.. Instrument Design and Performance for Optical/Infrared Ground-based Telescopes,
M. Iye, A. F. M. Moorwood, eds. (2003), vol. 4841 of Society of Photo-Optical Instrumentation Engineers
(SPIE) Conference Series, pp. 1657-1669.
63. J. A. Newman, et al., ApJS 208, 5 (2013).
64. M. C. Cooper, J. A. Newman, M. Davis, D. P. Finkbeiner, B. F. Gerke, spec2d: DEEP2 DEIMOS Spectral
Pipeline, Astrophysics Source Code Library (2012).
65. D. Lang, D. W. Hogg, D. J. Schlegel, ArXiv e-prints (2014).
66. A. E. Shapley, C. C. Steidel, M. Pettini, K. L. Adelberger, ApJ 588, 65 (2003).
67. R. Antonucci, ARAA 31, 473 (1993).
68. C. M. Urry, P. Padovani, PASP 107, 803 (1995).
69. D. B. Sanders, et al.,ApJ 325, 74 (1988).
70. P. J. McCarthy, ARAA 31, 639 (1993).
71. M. Lacy, et al., ApJS 154, 166 (2004).
72. D. Stern, et al., ApJ 631, 163 (2005).
73. W. N. Brandt, G. Hasinger, ARAA 43, 827 (2005).
74. N. L. Zakamska, et al., AJ 126, 2125 (2003).
75. M. Brusa, et al., ApJ 716, 348 (2010).
76. L. Hao, et al., AJ 129, 1783 (2005).
77. M. Kim, et al., ApJL 768, L9 (2013).
78. D. Stern, et al., ApJ 568, 71 (2002).
79. B. A. Groves, M. A. Dopita, R. S. Sutherland, ApJS 153, 75 (2004).
80. C. De Breuck, H. Rottgering, G. Miley, W. van Breugel, P. Best, A & A 362, 519 (2000).
81. T. Nagao, R. Maiolino, A. Marconi, A & A 447, 863 (2006).
82. K. N. Hainline, A. E. Shapley, J. E. Greene, C. C. Steidel, ApJ 733, 31 (2011).
83. M. Lacy, et al., ApJS 208, 24 (2013).
84. D. E. Vanden Berk, et al., AJ 122, 549 (2001).
85. T. Nagao, A. Marconi, R. Maiolino, A & A 447, 157 (2006).
86. W. Zheng, G. A. Kriss, R. C. Telfer, J. P. Grimes, A. F. Davidsen, ApJ 475, 469 (1997).
87. D. Stern, et al., ApJ 753, 30 (2012).
88. R. J. Assef, et al., ApJ 772, 26 (2013).
89. L. Yan, et al., AJ 145, 55 (2013).
90. M. Polletta, et al., ApJ 663, 81 (2007).
91. R. J. Assef, et al., ApJ 713, 970 (2010).
92. R. C. Hickox, et al.,ApJ 671, 1365 (2007).
93. R. C. Hickox, et al.,ApJ 731, 117 (2011).
94. Y. Shen, et al., AJ 133, 2222 (2007).
95. R. H. Becker, R. L. White, D. J. Helfand, ApJ 450, 559 (1995).
96. F. M. Montenegro-Montes, et al., MNRAS 388, 1853 (2008).
97. M. A. DiPompeo, M. S. Brotherton, C. De Breuck, S. Laurent-Muehleisen, ApJ 743, 71 (2011).
98. S. Fine, M. J. Jarvis, T. Mauch, MNRAS 412, 213 (2011).
99. P. J. E. Peebles, The large-scale structure of the universe (Princeton University Press, 1980).
100. J. N. Fry, ApJ 279, 499 (1984).
101. A. Cooray, R. Sheth, Physics Reports 372, 1 (2002).
102. P. F. Hopkins, G. T. Richards, L. Hernquist, ApJ 654, 731 (2007).
103. R. Reyes, et al., AJ 136, 2373 (2008).
104. E. Lusso, et al., ApJ 111 , 86 (2013).
105. E. Bertin, S. Arnouts, Astronomy and Astrophysics Supplement 117, 393 (1996).
106. D. B. Nestor, A. E. Shapley, C. C. Steidel, B. Siana, ApJ 736, 18 (2011).
107. R. A. C. Croft, G. B. Dalton, G. Efstathiou, W. J. Sutherland, S. J. Maddox, MNRAS 291, 305 (1997).
108. Y. Shen, et al., ApJ 719, 1693 (2010).
109. L. Guaita, et al.ApJ 733, 114 (2011).
110. K. L. Adelberger, C. C. Steidel, ApJL 627, LI (2005).
111. G. T. Richards, et al., AJ 124, 1 (2002).
112. R. van Ojik, et al., A & A 313, 25 (1996).
113. M. Villar-Martln, et al„ MNRAS 346, 273 (2003).
114. M. Villar-Martin, New A Rev. 51, 194 (2007).
115. A. Humphrey, et al., MNRAS 375, 705 (2007).
116. M. K. M. Prescott, C. L. Martin, A. Dey, ApJ 799, 62 (2015).
117. T. M. Heckman, G. K. Miley, M. D. Lehnert, W. van Breugel, ApJ 370, 78 (1991).
118. R. Trainor, C. C. Steidel, ApJL 775, L3 (2013).
119. G. Miley, C. De Breuck, A<L4 Rev. 15, 67 (2008).
120. K. K. Nilsson, J. P. U. Fynbo, P. Mollcr, J. Sommer-Larsen, C. Ledoux, A & A 452, L23 (2006).
121. D. J. B. Smith, M. J. Jarvis, MNRAS 378, L49 (2007).
122. T. M. A. Webb, et al„ ApJ 692, 1561 (2009).
123. J. W. Colbert, et al, ApJ 728, 59 (2011).
124. J. X. Prochaska, H.-W. Chen, A. M. Wolfe, M. Dessauges-Zavadsky, J. S. Bloom, ApJ 672, 59 (2008).
125. G. E. Prochter, J. X. Prochaska, J. M. O’Meara, S. Buries, R. A. Bernstein, ApJ 708, 1221 (2010).
126. G. J. Ferland, et al, PASP 110, 761 (1998).
127. A. Gould, D. H. Weinberg, ApJ 468, 462 (1996).
128. E. Lusso, et al, ArXiv e-prints (2015).
129. F. Arrigoni Battaia, et al., ArXiv e-prints (2014).
130. I. V. Strateva, W. N. Brandt, D. P. Schneider, D. G. Vanden Berk, C. Vignali, AJ 130, 387 (2005).
131. G. T. Richards, et al., ApJS 166, 470 (2006).
132. A. Baskin, A. Laor, J. Stern, MNRAS 438, 604 (2014).
133. F. Haardt, P. Madau, ApJ 746, 125 (2012).